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Abstract 

A new numerical discretization method for solving conservation laws is being devel- 
oped. This new approach differs substantially in both concept and methodology from 
the well-established methods-i.e., finite difference, finite volume, finite element, and spec- 
tral methods. It is motivated by several important physical/numerical considerations and 
designed to avoid several key limitations of the above traditional methods. 

As a result of the above considerations, a set of key principles for the design of numer- 
ical schemes was put forth in a previous report. These principles were used to construct 
several numerical schemes that model a 1-D time-dependent convection-diffusion equation. 
These schemes were then extended to solve the time-dependent Euler and Navier-Stokes 
equations of a perfect gas. It was shown that the above schemes compared favorably with 
the traditional schemes in simplicity ; generality, and accuracy. 

In this report, the 2-D versions of the above schemes, except the Navier-Stokes solver, 
are constructed using the same set of design principles. Their constructions are simplified 
greatly by the use of a nontraditional space-time mesh. Its use results in the simplest 
stencil possible, i.e., a tetrahedron in a 3-D space- time with a vertex at the upper time 
level and other three at the lower time level. Because of the similarity in their design, each 
of the present 2-D solvers virtually shares with its 1-D counterpart the same fundamental 
characteristics. Moreover, it will be shown that the present Euler solver is capable of 
generating highly accurate solutions for a famous 2-D shock reflection problem. Specifically, 
both the incident and the reflected shocks can be resolved by a single data point without 
the presence of numerical oscillations near the discontinuity. 



1. Introduction 


A new numerical discretization method for solving conservation laws is being devel- 
oped [1-8]. This new approach differs substantially in both concept and methodology 
from the well-established methods-i.e., finite difference, finite volume, finite element, and 
spectral methods [9-13]. It is conceptually simple and designed to overcome several key 
limitations of the above traditional methods. 

A two-level explicit scheme for solving a 1-D convection-diffusion equation was con- 
structed in [1] using this new method. Because the convection speed and the viscosity 
coefficient in the above equation are denoted by a and fi, respectively, the new scheme is 
referred to as the 1-D a-fi scheme. In [1], this scheme is subjected to a thorough theoretical 
and numerical analysis on stability, dissipation, dispersion, consistency, truncation error, 
and accuracy. The 1-D a-fi scheme is a two-way marching scheme [5], i.e., the forward 
marching scheme can be inverted to become the backward marching scheme. In other 
words, the marching variables at a lower time level can also be determined in terms of 
those at higher time levels. The special case of the the 1-D a-fi scheme which solves only 
the convection equation is referred to as the 1-D a scheme. This scheme is the only two- 
level explicit scheme known to the authors to be neutrally stable, i.e., free from numerical 
diffusion. Note that the solution of a pure convection problem has three fundamental prop- 
erties: (i) it does not dissipate with time, (ii) its value at a spatial point at a later time has 
a finite domain of dependence at an earlier time, and (iii) it is completely determined by 
the initial data at a given time. Because the a scheme is a two-level, explicit, and neutrally 
stable scheme, the above three properties are also shared by a solution of the a scheme. 
Contrarily, (i) a solution of a diffusive scheme will dissipate with time; (ii) the value of 
a solution of an implicit scheme at any space-time mesh point is dependent on all initial 
data, and all the boundary data up to the time level of the mesh point under considera- 
tion; and (iii) the unique determination of a solution of a multi-level scheme requires the 
specification of the initial data at two or more time levels. 

Because the a scheme is free from numerical diffusion and it is a special case of the 
a-p scheme with fi = 0, the a-p scheme has a special property that a classical scheme 
generally lacks, i.e., as the physical diffusion approaches zero, so does the numerical dif- 
fusion. Without this property, numerical dissipation may overwhelm physical dissipation 
and cause a complete distortion of solutions for problems with small viscosity. 

The 1-D a-fi scheme was derived again in [5] using a different type of conservation 
element and solution element. In [5], it also was extended to solve the 1-D time- dependent 
Navier-Stokes equations of a perfect gas. In spite of the fact that it does not use (i) 
any techniques related to the high-resolution upwind methods, (ii) any mesh-refinement 
techniques, (iii) any moving meshes, and (iv) any ad hoc parameter, the new Navier- 
Stokes solver is capable of generating highly accurate shock tube solutions. Particularly, 
for high-Renolds-number flows, shock discontinuities can be resolved almost within one 
mesh interval. 

The 1-D a scheme is neutrally stable and reversible in time. It is well known that such 
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a scheme generally becomes unstable when it is extended to model the Euler equations. It 
is also obvious that a scheme that is reversible in time c ann ot model a physical problem 
that is irreversible in time, e.g., an inviscid flow problem involving shocks. As a result, 
the 1-D a scheme was modified in [5]. Stability of this modified scheme, referred to as 
the 1-D a-e scheme, is limited by the CFL condition and 0 < e < 1 where e is a special 
parameter that controls numerical diffusion. The 1-D a-e scheme is reduced to the 1-D a 
scheme (which is free of numerical diffusion) when e = 0. 

The 1-D a-e scheme also was extended in [5] to become an Euler solver. The stability 
conditions of the Euler solver are similar to those of the 1-D a-e scheme. It was shown in 
[5] that the Euler solver is capable of generating accurate shock tube solutions within a 

wide range of CFL number. ' 

In this report, it will be explained how the 2-D versions of the a- pi, the a, and the 
a-e schemes can be constructed using the same set of design principles which were used to 
construct their 1-D counterparts. Because of this similarity in design, each of these 2-D 
versions virtually shares with its 1-D counterpart the same fundamental characteristics. 
Moreover, the 2-D a-e scheme will also be extended to become an Euler solver. This ex- 
tension is also very similar to its 1-D counterpart. We will not, however, discuss in this 
report the Navier- Stokes extension of the 2-D a- pi scheme. Because a Navier- Stokes prob- 
lem is fundamentally an initial-value/boundary value problem, the above extension (which 
is explicit and thus can not transmit information from one end of the boundary to another 
end in one time step) obviously cannot model such a problem unless the contribution of 
the viscous terms is small compared to that of the convection (inertial) terms. In general, 
this implies that the Navier-Stoke extension of the 2-D a- pi scheme is applicable only to 
high-Reynolds-number flows. In a future report, this extension will be introduced as a 
special case of a more general Navier-Stokes solver. 

The current method emphasizes simplicity, generality, and accuracy. It represents a 
clear break from the traditional methods in the basic concept of discretization. Most of 
the considerations that motivate its development and the key differences that separate the 
current method from the traditional methods were discussed in [5]. In the following, we 
present a more up -to -date version of these disc ussions: , 

(a) A set of physical conservation laws is a collection of statements of flux conservation in 
space-time. Mathematically, these laws are represented by a set of integral equations. 
The differential form of these laws is obtained from the integral form with the assump- 
tion that the physical solution is smooth. For a physical solution in a region of rapid 
change (e.g., a boundary layer), this smoothness assumption is difficult to realize by 
a numerical approximation that can use only a limited number of discrete variables. 
This difficulty becomes even worse in the presence of discontinuities (e.g., shocks). 
Thus, a method designed to obtain numerical solutions to the differential form with- 
out enforcing flux conservation is at a fundamental disadvantage in modeling physical 
phenomena with high-gradient regions. Particularly, it may not be used to solve flow 
problems involving shocks. Contrarily, a numerical solution obtained from a method 
that also enforces flux-conservation locally (i.e., down to a computational cell) and 
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globally (i.e., over the entire computational domain) will always retain the basic phys- 
ical reality of flux conservation even in a region involving discontinuities. For this 
reason, the enforcement of both local and globed £nx conservation in space and time 
is a tenet in the current development. The concept of space- time conservation element 
is introduced to serve this purpose. 

Among the traditional methods, finite difference, finite element, and spectral 
methods are designed to solve the differrential form of the conservation laws. Note 
that the set of integral equations usually solved in a finite-element scheme is equiv- 
alent to the differential form of the conservation laws assuming certain smoothness 
conditions. However, these integral equations generally are different from the integral 
equations representing the conservation laws. Even if they are cast into a conservative 
form, the resulting flux- conservation conditions generally do not represent the physical 
conservation laws. 

The finite volume method is the only traditional method designed to enforce flux 
conservation. A finite- volume scheme may enforce flux conservation in space only, 
or in both space and time. As a preliminary to this enforcement, a flux must be 
assigned at any interface separating two neighboring conservation cells. In a typical 
finite- volume scheme, it is evaluated by extrapolating or interpolating the mesh values 
at the neighboring cells. This evaluation generally requires an ad hoc choice of a 
special flux model among many models available [14-16]. Generally numerical results 
obtained are dependent on which model one chooses. Also this process of interpolation 
and extrapolation generally is time consuming and may result in numerical smearing. 

Contrarily, by using the concept of space-time solution element, and considering 
the spatial derivatives of dynamic variables as independent variables (to be discussed 
further shortly), in the current method, flux evaluation at an interface is carried out 
without interpolation or extrapolation. It is an integral part of the solution procedure, 
(b) The numerical variables used in a spectral method, i.e., the expansion coefficients, 
are global parameters pertaining to the entire computational domain. As a result, a 
spectral method generally (i) lacks local flexibility and thus may be applied only to 
problems with simple geometry, and (ii) is hindered by the fact that it must deal with 
a full matrix that is difficult to invert. 

By design, only local parameters will be used in the current method. Moreover, the 
set of discrete variables in any one of the numerical equations to be solved generally is 
associated with only two neighboring solution elements. The exception to this general 
rule occurs only in the situation in which numerical diffusion is to be introduced 
deliberately (see [5], [6], and Secs. 3 and 4). Even in this special case, only the discrete 
variables associated with a few immediately neighboring solution elements will enter 
any equation to be solved. Thus, one needs only to deal with a very sparse matrix. 
Moreover, the maximum number of solution elements involved in a numerical equation 
of the current discretization framework is independent of the order of accuracy of a 
particular scheme. Contrarily, the order of accuracy of a classical finite-difference 
scheme generally can be increased only by using variables at more mesh points in 
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each of its equations. Usually, a side effect of this practice is an increase in numerical 
diffusion. Also it may be difficult to implement a high-order finite-difference scheme 
near a boundary because there are no real mesh points outside the boundary. Note 
that, in the absence of body force, direct physical interactions occur only among the 
immediate neighbors. The current design is also consistent with this physical reality. 

(c) Space and time traditionally axe treated separately in the time marching schemes. 
Generally one obtains a system of ordinary differential equations with time being the 
independent variable after a spatial discretization. As an example, elements in the 
finite element method usually are used for spatial discretization. These elements axe 
domains in space only. 

Because flux conservation is fundamentally a property in space-time, space and time 
axe unified and treated on the same footing in the current method. Thus, conservation 
elements and solution elements used in the time-dependent version of the current 
method are domains in space-time. The significance of this unified approach cannot 
be overemphasized. It makes it easier for a numerical analogue to share the same 
space-time symmetry of the physical laws (see [2], [5], [6] and the following sections). 

(d) In a finite-difference scheme, derivatives at mesh points are expressed in terms of mesh 
values of dependent variables by using finite-difference approximations. Accuracy of 
these approximations, especially those of higher-order accuracy, generally is excellent 
as long as dependent variables vary slowly across a mesh interval. However, it may 
not be adequate if these variables vary too rapidly. Thus, in a high-gradient region, 
e.g., a boundary layer, accuracy may demand the use of an extremely fine mesh. In 
turn, a prohibitively high computing cost may result. 

The current method avoids the above pitfall by expressing the numerical solu- 
tion within a solution element as an expansion in terms of certain base functions. 
As in a spectral method, the expansion coefficients are considered as the indepen- 
dent numerical variables to be solved for simultaneously. For simplicity, Taylor’s 
expansions will be used in the current paper. For this special case, the expansion 
coefficients are interpreted as the numerical analogues of the derivatives. Note that: 

(i) Because the derivatives are considered as independent variables, their values at 
the initial/boundary surfaces may be specified as a part of the initial /boundary con- 
ditions. This specification may result in more accurate initial/boimdary conditions; 

(ii) van Leer [17] also has attempted to improve accuracy by introducing two indepen- 
dent numerical variables for each independent physical variable, and (iii) the c ur rent 
solution procedure has no resemblance with those used in compact difference schemes. 

(e) With a few exceptions, numerical diffusion generally appears in a numerical solution 
of a time-marching problem. In other words, the numerical solution dissipates faster 
than the corresponding physical solution. For a nearly inviscid problem, e.g., flow with 
a high Reynolds number, this could be very serious because numerical dissipation may 
overwhelm physical dissipation and cause a complete distortion of solutions. One may 
argue that numerical diffusion can be reduced by increasing the order of accuracy of 
the scheme used. However, because the order of accuracy of a scheme is generally 
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determined with the aid of Taylor’s expansion, and the latter is valid only for a 
smooth solution, it has meaning only for a smooth solution. Thus the use of a scheme 
of higher-order accuracy may not reduce numerical diffusion associated with high- 
frequency Fourier components of a numerical solution. This is the reason that the 
Leapfrog scheme, which is free from numerical diffusion, can outperform schemes with 
higher-order accuracy in solving some wave equations [18]. 

In a study of finite-difference analogues of a simple convection equation [2], it was 
shown that a numerical analogue will be free from numerical diffusion if it does not 
violate certain space-time invariant properties of the convection equation. In other 
words, numerical diffusion may be considered as a result of symmetry-breaking by 
the numerical scheme. Because of its intrinsic nature of space-time unity, the current 
framework is an excellent vehicle for constructing a numerical analogue that shares 
the same space- time invariant properties with the physical equation (see [5], [6], and 
the following sections). 

It is recognized that a certain amount of numerical diffusion may be needed to pre- 
vent large dispersive errors [19] that are often caused by the presence of high-frequency 
disturbances (such as round-off errors). Therefore, schemes were constructed such that 
the numerical diffusion can be controlled by a single adjustable parameter (see [5], [6], 
and Secs. 3 and 4). The numerical diffusion is shut off when this parameter is set to 
zero. 

(f) High- resolution upwind methods [13], which we consider to be a branch of the finite 
volume method, are heavily dependent on characteristics-based techniques. For the 
1-D time-dependent case, the characteristics are curves in space-time, and the coeffi- 
cient matrix associated with the Euler equations [20] also can be diagonalized easily. 
As a result, these techniques are easy to apply. However, for multidimensional cases, 
the characteristics are 2-D or 3-D surfaces in space-time [21]. Moreover, the coefficient 
matrices cannot be diagonalized simultaneously by the same matrix [20]. Because of 
the above complexites, application of these techniques to multi-dimensional problems 
is much more diffucult. Furthermore, high-resolution methods generally require the 
use of ad hoc parameters, e.g., flux-limiters and/or slope-limiters, and other ad hoc 
techniques. These ad hoc techniques may lead to numerical diffusion which varies 
from one place to another and from one Fourier component to another. In other 
words, numerical solutions may suffer annihilation of sharply different degrees at dif- 
ferent locations and different frequencies [22]. Also, these techniques generally are 
also difficult to extend to a space of higher dimension. 

Because the current framework is developed to solve multidimensional problems 
(see the following sections), simplicity and generality weigh heavily in its design. Thus, 
we do not use the characteristics-based techniques, and also try to avoid the use of ad 
hoc techniques. Moreover, the concept of characteristics generally is not applicable 
to the Navier-Stokes equations, which are non-hyperbolic in nature. Therefore, the 
above decision also makes it easier for the current framework to solve the Navier-Stokes 
equations [5]. 
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This completes the discussion of the motivation for the current development. In sum- 
mary, the development is guided by the following requirements: (i) To enforce both local 
and global flux consevation in space and time with flux evaluation at an interface being 
an integral part of the solution procedure and requiring no interpolation or extrapolation; 
(ii) To use local discrete variables such that the set of variables in any one of the numer- 
ical equations to be solved is associated with a set of immediately neighboring cells; (iii) 
Space and time are unified and treated on the same footing; (iv) Mesh values of dependent 
variables and their derivatives are considered as indep enden t variables to be solved for 
simultaneously; (v) To minimize numerical diffusion, a numerical analogue should be con- 
structed, as much as possible, to be compatible with the space-time invariant properties 
of the corresponding physical equations; and (vi) To exclude the use of the characteristics- 
based techniques, and to avoid the use of ad hoc techniques as much as possible. It is 
the purpose of this report to show that the above requirements can be met with a simple 
unified numerical framework even for multi dim ensional problems. 
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2. The a-fi Scheme 


In this section, we consider a dimensionless form of the 2-D convection-diffusion equa- 
tion, i.e., 


du du du ( 6Pu dPu\ 

Tk +a ’7h* a, fy-^\J^ + lwj =a 


( 2 . 1 ) 


where a x , a y , and p (> 0) are constants. Let x x = x, x 2 — y, and x 2 = the the coordinates 
of a three-dimensional Euclidean space E 3 . By using Gauss’ divergence theorem in the 
space-time E 3 , it can be shown that Eq. (2.1) is the differential form of the integral 
conservation law 

(f h • ds — 0. (2.2) 

Js(V) 

Here (i) S(V ) is the boimdary of an arbitrary space-time region V in E3, (ii) 


r def 


h = ( a x u — pdu/dx, a y u — pdu/dy , u) 


(2.3) 


is a current density vector in E 3 , and (iii) ds — daft with da and n, respectively, being the 
area and the outward unit normal of a surface element on S'(V’). Note that (i) h ■ ds is the 
space-time flux of h leaving the region V through the surface element ds, and (ii) all math- 
ematical operations can be carried out as though E 3 were an ordinary three-dimensional 
Euclidean space. As will be shown shortly, E 3 will be divided into nonoverlapping space- 
time regions referred to as conservation elements (CEs). 

Let n denote the time level and 

f nd =nAt, n = 0, ±1/2, ±1, ±3/2,.... (2.4) 


Let j and k be spatial mesh indices with j, k = 0, ±1/3, ±2/3, ±1,... (see Figs. 1- 
4). Let fli denote the set of mesh points (j,k,n) with j,k = 0, ±1,±2, ..., and n = 

±1/2, ±3/2, ±5/2, These mesh points are marked by solid circles. Let denote the 

set of mesh points (j, k, n ) with j, k = 1/3, 1/3 ± 1, 1/3 ± 2, . . ., and n = 0, ±1, ±2, — 
These mesh points are marked by open circles. The union of S7i and 1^2 will be denoted 
by Q. 

Each mesh point (/, k, n+1/2) in £2i (by definition, this implies that n = 0, ±1, ±2, . . .) 
is associated with three CEs, denoted by CE^(j, k, n + 1/2), £ = 1,2,3 (see Fig. 5(a)). It is 
also associated with a solution element (SE), denoted by SE^^(jf, k, n + 1/2) (see Fig. 5(b)). 
Similarly, each mesh point (j, k, n + 1) in 0 2 is associated with three conservation elements 

CE ^\j,k,n + 1), £ — 1,2,3 (see Fig. 6(a)), and a solution element SE ( 2 \j,k,n + 1) (see 
Fig. 6(b)). Each CE is a quadrilateral cylinder in space- time while each SE is the union of 
three vertical planes, a horizontal plane, and their immediate neighborhood. The geometry 
of the hexagon ABCDEF , which appears in both Figs. 5(a) and 6(a), is determined by 
three positive parameters w, b and h (see Fig. 7(a)). Without any loss of generality, we 
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assume that the line segment joining points D and A in Fig. 7(a) is parallel to the x-axis. 
Note that the form of Eq. (2.1) will not change under an orthogonal transformation on the 
x-y plane. Thus one always can introduce a set of Cartesian coordinates (x,y) such that 
the above line segment is parallel to the x-axis. However, the values of a x and a y may 
change as a result of such a transformation. 

According to Fig. 1, E 3 can be filled with the CEs defined above. Moreover, it is seen 
from Figs. 5(a), 5(b), 6(a), and 6(b) that the boundary of a CE is formed by the subsets 
of two neighboring SEs. 

Let the space- time mesh be uniform, i.e., the parameters At, w, b, and h be constants. 
Let Xj t fc and j/j,* be the x- and y- coordinates of any mesh points (j, k, n) € fl. Let xo,o = 0 
and yo t o = 0. Then information provided by Figs. 7(a) and 7(b) implies that 


x i,k = O' + k)to + (k- j)b, y jtk = (k - j)h. (2.5) 

Let Hi, n 2 , n 3 , n 4 , n 5 , and n 6 be the vectors depicted in Fig. 7(a). They lie on the x-y 
plane and are the outward unit normals to AB, BC , CD, DE , EF, and FA, respectively. 
It can be shown that 



(h,-b+w/ 3,0) 


(2.6a) 

Hi = 

v // l 2 + (6_ u ,/3)2’ 

n 4 = —Hi, 


n 2 = (0,1,0), 

pi 

Cn 

II 

1 

5* 

(2.66) 

n 3 = 

(— ft, 6-f ic/3, 0) 

VV + (6 + u>/3) 2 ’ 

«6 = -n 8 , 

(2.6c) 


For any (j, k, n) £ 0, let 


SE(j, k, n) = 


SE (1) (j,fc,n), if ( j,k,n ) € 
SE (2) (j, k,n), if (j, k,n) 6 fi 2 . 


(2.7) 


For any (x,y,t) € SE {j, k,n ), u(x, y,t) and h(x,y,t), respectively, are approximated by 


M) = «>,* + («*);,*(* - x i>*) + C «>)?,*(» - vm) + («*)",*(* - <n )> ( 2 - 8 ) 


and 


h*(x,y,t-,j,k,n) d = [a r u*(x,y,<;j,fc,n) - fxdu*(x,y,t-,j,k,n)/dx, 
y, t'li, k, n ) - fidu*(x, y, t ; j, k, n)/dy, u*(x, y,<; fc, n)] , 


(2.9) 


where u” t , (u*)” fc , (u y )" fc , and (itt)£ fc are constants within SE(j, fc, n). The last four 
coefficients, respectively, can be considered as the numerical analogues of the values of u, 


8 


iliiliilllliil 



du/dx, du/dy , and du/dt at ( Xj,yk,t n ). As a result, the expression on the right side of 
Eq. (2.8) can be considered as the first-order Taylor’s expansion of u(x,y,t ) at ( Xj,yk,t n )• 
Also note that Eq. (2.9) is the numerical analogue of Eq. (2.3). 

We shall require that u = u*(x,y,t-,j,k,n ) satisfies Eq. (2.1) within SE(j, k, n). As a 
result, 

( u t)j,k ~ ~ [ a *( w r)£* + fly( u ff)j,Jt] • (2.10) 

Because Eq. (2.8) is a first-order Taylor’s expansion, the diffusion term in Eq. (2.1) has 
no counterpart in Eq. (2.10). As a result, the diffusion term has no impact on how 
u*(x,y,t]j,k,n ) varies with time within SE(j,k,n). However, as will be shown shortly, 
through its role in the numerical analogue of Eq. (2.2), it does influence time- dependence 
of numerical solutions. Note that, for a higher-order scheme, how u*(x, y,t;j, k,n ) varies 
with time within SE (j, k, n) will be influenced by the presence of the diffusion term. Sub- 
stituting Eq. (2.10) into Eq. (2.8), one has 

u*(x,y,t;j, k y n) = + («*)"* [(* - x jtk ) - a x (t - t n )] 

+ K)",* Kv - y;,*) - <*y(t - <")] • 


Thus there are three independent marching variables, i.e., (u r )"fc> and (uj ,)" k asso- 
ciated with a mesh point ( j,k,n ) £ Cl. For any ( j , fc,n + 1/2) £ ftj, these variables will 
be determined in terms of those associated with the mesh points ( j + 1/3, k + 1/3, n), 
( j — 2/3, k + 1/3, n), and (J + 1/3, fc — 2/3, n) (see Fig. 8(a)) by using the flux conservation 
relations: 

(f h*-ds = 0, l — 1,2,3. (2.12) 

JS(CE?\j,k,n+ 1 / 2 )) 

Similarly, the marching variables at any ( j , k, n + 1) £ O 2 are determined in terms of those 
associated with the mesh points ( j — 1/3, Jfc + 2/3, n + 1/2), (j — 1/3, k — 1/3, n + 1/2), and 
(j + 2/3, k — 1/3, n + 1/2) (see Fig. 8(b)) by using the flux conservation relations: 


/ h*-ds = 0, £= 1,2,3 . 

JS(CEf\j,k,n+ 1 )) 


(2.13) 


Obviously, Eqs. (2.12) and (2.13) are the numerical analogues of Eq. (2.2). 

As a result of Eqs. (2.12) and (2.13), the total flux leaving the boundary of any CE 
is zero. Because the flux at any interface separating two neighboring CEs is calculated 
using the information from a single SE, the flux entering one of these CEs is equal to that 
leaving another. It follows that the local conservation conditions Eqs. (2.12) and (2.13) 
will lead to a global conservation condition, i.e., the total flux leaving the boundary of any 
space-time region that is the union of any combination of CEs will also vanish. 

In the following, several preliminaries will be given prior to the evaluation of Eqs. (2.12) 
and (2.13). To proceed, note that a mesh fine with j and n being constant or a mesh line 
with k and n being constant is not aligned with the z-axis or the y-axis. We shall introduce 
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a new spatial coordinate system (£, t/) with its axes aligned with the above mesh lines (see 
Fig. 7(c)). 

Let e x and e y be the unit vectors in the x- and the y- directions, respectively. Let 

e*£ and e v be the unit vectors in the directions of Djfr and DB (i.e., the j- and the k- 
directions-see Figs. 7(a)-(c)), respectively. It can be shown that 


and 

where 

and 


e c = [(w - b)e x - ht y ) /a(, 
e, = [(u> + b)e x + he y ] /at/, 

aC =' |5?| = - i) 1 + h 2 , 

A17 d = \DB\ = y/(w~+ 6) 2 + ft 2 . 


(2.14) 

(2.15) 

(2.16) 
(2.17) 


Let the origin of (x, y) also be that of (£, rj). Then the spatial coordinates (x, y) and ((, J7) 
of any point in E$ are related by the condition 


Ce c + r?e n = ie, + ye„. 

Substituting Eqs. (2.14) and (2.15) into Eq. (2.18), one has 


(;)-'©■ 


(2.18) 


(2.19) 


( 2 . 20 ) 


Here 


def 


T~1 def 


' w — b 

w + b\ 

aC 

AT/ 

h 

h 

. a< 

AT/ / 

aC 

(w + 6 )aC 

2w 

2 wh 

A»7 

(w — b) AT/ 

2 w 

2wh 


( 2 . 21 ) 
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Note that the existence of T 1 , the inverse of T, is assured if wh ^ 0. 

With the aid of Eqs. (2.5), (2.20), and (2.22), it can be shown that the coordinates 
((, 77 ) of any mesh point (j, k,n) E Cl are given by 

C = j a £, and r] = k at), (2.23) 

i.e., a( and A rj are the mesh intervals in the £- and the rj- directions, respectively. 

Next we shall introduce several coefficients which are tied to the coordinate system 


(C,rj). Let 

d JL f T -l 

(:;)■ 

(2.24) 

Also, for any ( j , k,n) € Cl, let 

I 



(2.25) 

\ 

K( u v)lk) 




where T* is the transpose of T. For those who are familiar with tensor analysis, the 
following comments will clarify the meaning of the above definitions: 

(a) (a^ ,a v ) are the contravariant components with respect to the coordinates ((, rj) for 
the spatial vector whose x- and y- components are a x and a y , respectively. 

(b) ((«{)£*, (ii^)” fc ) are the covariant components with respect to the coordinates ((,tj) 
for the spatial vector whose x- and y- components are (u r )" t and (%)"*, respectively. 

(c) Because the contraction of the contravariant components of a vector and the covariant 
components of another is a scalar, Eq. (2.10) can be rewritten as 

(«<)",* = - [ a <( w <)" ( * + a v( u v)],k] • (2.26) 


(d) Under the linear coordinate transformation defined by Eqs. (2.19) and (2.20), (C — 
jA^,rj — k at]) are the contravariant components with respect to the coordinates ((, tj) 
for the spatial vector whose x- and y- components are x—Xj t * and y—yj,k, respectively. 
Using the same reason given in (c), Eq. (2.11) implies that 

u*(x,y,t;j,k,n) = u*((,T),t-,j,k,n), (2.27) 


where 


«*(C» *; h n) d = + («<)”* [(C - j a C) - <*<(* ~ *")] 

+ ( u v)j,k l(»7 - k *v) ~ <**(* - t n )] . 


(2.28) 


Note that Eqs. (2.26) and (2.27) can also be verified directly using Eqs. (2.20), (2.22), 
(2.24), and (2.25). 
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Next, let (i) 



j. def 6 

and 

4. def 6 

s = —<V> 

1 at? 

(2.29) 

(ii) 


and 

v u v )>>* g y u v)j,ki 

(2.30) 

(iii) 

def At 4. 
v < = 

and 

def At i 

v * = T a ^‘ 

(2.31) 

and (iv) 





* def 

u = 

9flAt / >\2 t def 9/iAt 
8w 2 h 2[ W ^ " 8 iv 2 h 2 

-(*r?) 2 , 

, . def 9 /i At . , 2 

and = 8to 2 /i 2 AT ’ 

(2.32) 

where 

def 
AT = 





2 V& + h 2 . 

(2.33) 


The coefficients defined in Eqs. (2.29) and (2.30) can be considered as the normalized 
counterparts of those defined in Eqs. (2.24) and (2.25). Also note that a(, at?, and at, 
respectively, axe the lengths of the three sides DF, BD , and FB of A BDF depicted in 



Figs. 7(a)-(c). Moreover, by substituting Eq. (2.29) into Eq. (2.31), one has 




at At 2 , a„ At 2 

aC = r< , and a>) =-n,. 

(2.34) 

= 


In other words, (2/3)^ and (2/3)u r) are the Courant numbers in the £- and 77 - 
respectively. 

directions, 



Furthermore, let cr^ - ,. . ., be defined by 




( 1 )± def - 
°11 = 1 - 

(2.35) 

I 


^S }± d = it 1 - - v v){l + v c) + in + ir - £c> 

(2.36) 


- 

* 1 ?* djf _ „ c _ „ f )(i + „ w ) + + £ r - 

(2.37) 

1 ? 

i 

( 1 )± def . . 

^21 = 1 + V<1, 

(2.38) 



= T (1 + r <)(2 - n c ) - 2 £„ 

(2.39) 

m 


<*23^ *= ±(1 + I/ <)(1 + + Cc + (v £t, 

(2.40) 

■ 

a 

( 1 )± def - . 
*31 = 1 + 

(2.41) 


- 

12 


- 



^ 2 )=t = f ±(* + ^)(! + v d + U + - £r, 

(2.42) 

ff 33 ^ = “F(l + I/ I?)(2 — I'j;) — 2^, 

(2.43) 

( 2 )± def . 

*11 = 1 + z/ c + I/„ 

(2.44) 

off * d = =F(1 + t'c + - ^c) -tv-tr + U’ 

(2.45) 

off* drf + u< + Ufi ^i _ l/jj ) - £«. - £ r + 

(2.46) 

( 2 )± def - 
*21 = 1-^1 

(2.47) 

off* def + + 

(2.48) 

*ff * = T (1 - I/ c )( 1 - *,) - - 6 , + €r, 

(2.49) 

( 2 )± def 

O 31 = l-^i 

(2.50) 

*32^ = T (1 — I/ i?)(l — y c) — — Ci? + £r s 

(2.51) 

off* d = ±(1 - i /,)(2 + I /,) + 2 ^ C . 

(2.52) 


Note that: 

(a) Each of Eqs, (2.35)-(2.52) represents two equations. One corresponds to the upper 
signs while another, to the lower signs. 

(b) The definitions given in Eqs. (2.35)-(2.43) will be used in the first marching step of 

the a-jj, scheme; while those given in Eqs. (2.44)-(2.52) will be used in the second 
marching step. It is seen that the expressions on the right sides of the former can 
be converted to those of the latter, respectively, by reversing the “+” and ” signs. 
Moreover, for every pair of m and £, and are converted to and , 
respectively, if z/£, and are replaced by -i^, -v v , — and -f r , 

respectively. 

Equations (2.12) and (2.13) are evaluated in Appendix A. With the aid of the above 
definitions, the results are summarized as follows: 

(a) Eq. (2.12) with t— 1: 


o-riH , (i)+ + i o-ri)+ +1 n+1 / 2 

a n u -f- <t 12 t <t 13 u v 
(b) Eq. (2.12) with t = 2: 


^ + off 


Jj+l/3,Jfc+l/3 


(2.53) 


r ( D+ 

'21 


u + a. 


( i )+„+ a . / T ( i )+„+ r +1/2 _ 


22 Uj + < 7 . 


23 u n 


J,* 


-(1)- . _(!)-„+ , _(!)-„+' 

^21 U ' ^22 T °23 U »7 


j-2/3,fc+l/3 


(2.54) 
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(c) Eq. (2.12) with £ = 3: 


4V + « + 4? + «c + 4s )+ “?] > k = [ a 3i u + 4i u t + 4z u j] y+1/3 k _ 2/3 • 

3 ’ 3 (2.55) 

(d) Eq. (2.13) with £ = 1: 


'&)+„ 4. /t ( 2 >+, / + 4. _(2)+ + l n+1 _ (2)- 

a ll U + <T 12 + <T 13 U v — <T n U 


+ U C + a \3 


, n+1/2 

K 

11 ij-l/3,k- 


(e) Eq. (2.13) with £ = 2: 


1/3 

(2.56) 


L( 2 )+ u + < t (2)+ « + 4- (7 (2)+ ti + l n+1 - Lt (2)- u + < t (2 >"u+ 4- ^-u+1 " +1/2 
["21 u + "22 U C + °23 — [ a 21 u ^ a 22 ' °23 Uf ? J J+2 / 3 Jfc_ 


(£) Eq. (2.13) with £ = 3: 


1/3 

(2.57) 


[ cr 3i )+u + "?2 )+u c + ct 33 )+u J]^ = [4 2 1 « + 4? «c+ag } ttj] j _ i/ / 3fc+2/3 - 

J ’ (2.58) 

Here (j, k,n + 1/2) £ fix is assumed in Eqs. (2.53)-(2.55); while (j, n + 1) £ fl 2 is assumed 
in Eqs. (2.56)-(2.58). Also, to simplify notation, in the above and hereafter we adopt a 
convention that can be explained using the expression on the left side of Eq. (2.56) as an 
example, i.e., - L " • 

_(2)+ . _( 2 )4 4 . (2)4 4 ] " +1 def L(2)+ n+1 I /T (2)+/.,4\n41 i /T ( 2 ) + / 7# + ) n +ll 
"ll U + a 12 “c ' "l3 “ijj^ ~ [°il u j,k + "l2 l u < )j,k +°i3 \ U V )j,k J- 


Consider the special case with /i = 0. It follows from Eq. (2.32) that 

£< = 6,=£r=0. (/x = 0) 


(2.59) 


Combining Eq. (2.59) with Eqs. (2.35)-(2.37), one concludes that cr^*, and 

contain a common factor (1 — — v^). Similarly, each of three consecutive pairs of 

coefficients defined in Eqs. (2.38)-(2.52) also contain a common factor. As a result, one 
concludes that: 

(a) Eq. (2.53) is satisfied by either 1 — vq — 1 /, = 0 or 


■ n+l/2 , 1 n 

U + (1 + I/ C )«+ + (1 + !/,)«+ = U - (1 + I/ C )uJ - (1 + U V )U+ 

J j f k l J j-f-l/3, 


(b) Eq. (2.54) is satisfied by either 1 + = 0 or 


*41/3 

(2.60) 


r , -in41/2 r i ,i» 

u - (2 - I/ C )u+ + (i + i/,)u+ = u + (2 - i/<)uj - (1 + u v )u+\ 

L s j,k s J j — 2/3, 


*41/3 

(2.61) 
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(c) Eq. (2.55) is satisfied by either 1 + u v = 0 or 

| n+l/2 


r , -in+i/z , 

u + (1 + ^c)“< - (2 - .. = « - (1 + ^c) w c + ( 2 - 


(d) Eq. (2.56) is satisfied by either 1 + + i/, = 0 or 

in+l 


J>+l/3,*-2/3 

(2.62) 


r i 1 1 n+1 r , . n+1/2 

[« — (1 — I/ C )u+ - (1 - V n )u+\ ' k = [ti + (1 - V < )ttj + (1 - J i _ 1/3iJfc _ 1/3 ■ 

h 3 ’ (2.63) 

(e) Eq. (2.57) is satisfied by either 1 — = 0 or 


u + (2 + i^)u£ ~ C 1 ~ u v) u 


n+l 

i.fc 


u - (2 + i/ c )uj + (1 - v v )u+ 


(f) Eq. (2.58) is satisfied by either 1 — = 0 or 

in+l 


-] n+l/2 
j+2/3,*-l/3 ’ 

(2.64) 


n-t-i r , 1 n+l/2 

u-(l-z/ c )w+ + (2 + i/,)u+J^ = [ti + (l-i/ c )ttJ-(2 + i/,)«i+J^ i ^ t+2/3 

(2.65) 

Here (j, k , n+l/2) € is assumed in Eqs. (2.60)-(2.62); while (j, k, n + l) 6 is assumed 
in Eqs. (2.63)-(2.65). The current “a” scheme, i.e., the scheme that solves Eq. (2.1) with 
fi = 0, will be constructed using Eqs. (2.60)-(2.65) instead of Eqs. (2.53)-(2.58). Assuming 
H — 0, Eqs. (2.60)-(2.65) imply Eqs. (2.53)-(2.58). However, the reverse is false unless one 
assumes that 

[1 - (i/ c + u v ) 2 ] (1 - v\) (l - i/J) / 0. (2.66) 

Note that the expressions within the brackets in the first three equations in Eqs. (2.60)- 
(2.65), respectively, can be converted to those in the last three by reversing the “+” and 
— signs. 

Let sj^, s^\ and denote the expressions on the right sides of 

Eqs. (2.60)-(2.65), respectively. Then it can be shown that Eqs. (2.60)-(2.62) are equiva- 
lent to 

“"J 1/2 = | [f 1 - "C - + (1 + v d s< 2 > + (1 + <'»)»» 1) ] > 


-*?>). 


and 


where ( j,k,n + 1/2) € fli- Also Eqs. (2.63)-(2.65) are equivalent to 

“i.t 1 = \ [(! + + v v) s i ) + C 1 ~ + (1 - v v) s z 2) 


(2.67) 

( 2 . 68 ) 

(2.69) 

(2.70) 
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and 


5 


( 2 . 71 ) 
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Let the 3x3 matrices of rank one Q\ , k = 1,2, and £ — 1, 2, 3, be defined by 


, (2-80) 

where the row matrix (ji is the transpose of the column matrix Note that 

the explicit forms of the matrices defined above, respectively, can be obtained from 
Eqs. (3.50)-(3.55) by letting e = 0. 

By their definitions, each of s\ k \ k = 1,2, and £ = 1,2,3, can be expressed as the 

product of one of the row matrices , k = 1 , 2 and £ = 1, 2, 3, and one of the column 

matrices defined in Eq. (2.73). As an example, 

s[ 1] = ($ 1} ) ‘ q(j + 1/3, k + 1/3, n). (2.81) 

Thus, Eqs. (2.67)-(2.69) can be expressed as 


(2.81) 


q O', k, n + 1/2) = Q[ 1} q(j + 1/3, k + 1/3, n) + Q^q (j - 2/3, k + 1/3, n) 


(2.82) 


+ Q^q(j + 1/3, A: — 2/3, n), ( j,k,n + 1/2) G 

Similarly, Eqs. (2.70)-(2.72) can be expressed as 

q (i, k, n + 1) fe Q { ?q (j - 1/3, k - 1/3, n + 1/2) + Q™q(j + 2/3, k - 1/3, n + 1/2) 

+ Q?qU - 1/3 ,k + 2/3, n + 1/2), (j,k,n+ 1) G ^ 2 - 

(2.83) 

The marching procedure in the a scheme is formed by applying the marching steps defined 
by Eqs. (2.82) and (2.83) successively. 

As a preliminary for future development, we apply Eq. (2.82) and then Eq. (2.83). 
The result is: 

q(j,k,n + 1/2) = Q^’Q^’qij + l,*,n - 1/2) + Qj’ q(],k + l,n — 1/2) 

+ O' - 1 , k,n - 1/2) + Q^Qfqd - 1 , k + 1 , n - 1/2) 

+ U, * - l.» - 1/2) + Qi^Q^qU + 1, i - 1, n - 1/2) 

+ WiVi 2 ' + Q^Q? + Q?Q?)qU, k,n- 1/2), 

(2.84) 

where ( j,k,n + 1/2) G Oi. Similarly, by applying Eq. (2.83) and then Eq. (2.82), one 
concludes that 

q(j, k, n + 1) = Q^Q^qU - 1, k, n) + Q^Q^qij, k - 1, n) 

+ Q^Q^qU + 1, k, n) + Q^Qi^qU + l,k-l,n) 

+ Q ( 3 ) Q? ) q(h k + l,n) + Q^Q^q O' - 1, k + 1, n) 

+ + Q^Q^qU, M), 
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where ( j , k, n + 1) 6 H 2 .' 

Next we shall consider the flux conservation conditions Eqs. (2.53)-(2.58) assuming 
fj > 0. Let and be the determinants of the matrices formed by <rj^ + and <rj^ + , 
k,£ = 1,2,3, respectively. It is shown in Appendix A that, for any /i, 


A (1) = 3 


3(1 + i/ c )(l + i/,)(l - i/ c - !/„) + 2(1 + i/ c )(l -v<- i/„)f c 


+ 2(1 + !/,)(! - v c - + 2(1 + i/ c )(l + + 



( 2 . 86 ) 


and 


a< 2 > = 3 


3(1 - I^Xl - l/„)(l + v < + u„) + 2(1 - I/ C )(l + v< + Vr,)^ 


+ 2(1 — t'ij)(l +V£ + Vv)£ri + 2(1 — ^)(1 — V^r + 



(2.87) 


Because (i > 0, ^ > 0, > 0, and £ r > 0, it follows from Eqs. (2.86) and (2.87) that 

A^ > 0 and A*' 2 ) > 0 if 


M < 1, Wn\ <1, and \v{ + i/ n \ < 1. (2.88) 

In the following dicussions, we assume A^ / 0 and A^ ^ 0. 

Let i^, ij 1 ^, a£\ s^\ s£\ and denote the expressions on the right sides of 
Eqs. (2.53)-(2.58), respectively. Let 

Pn = [3(1 + ^<)(1 + + (3 + 2u^ + Vrj)U + (1 + v v)((v - £r)] /a (1 \ (2.89) 

P 22 = [-3(1 + V,)(l - V< - v n ) - (3 - 2i / c - *,)£ + (1 + */„)(£„ - e r )l /a (1) , (2.90) 

$ = [(^c + + (2 - *,X*r - e, j] /* (1) , (2.91) 

/4V =* [3(1 + i / c)(l + v , ) + (3 + + 2v n )£ ri + (1 + ^c)(^C — Cr)] /a (1) , (2.92) 

4V = [(*< + 2v v )i v + (2 - iO(* T - €c)J /A®, (2.93) 

4V d = [-3(1 + ^)(1 - *>< - ^) - (3 - v< - 2i/,X, + (1 + */ c X€c “ &•)] /A (1) , (2.94) 

4? d = [-3(1 — ^<)(1 — Vj,) - (3 - 2i/ c - i/„)£ c + (1 - ^)(^r - £,)] /a (2) , (2.95) 

4V d = [3(1 — ^i?)(l + ^ + (3 + 2u<; + i/,)^ + (1 — ^)(^r — &»)] /a (2 \ (2.96) 
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(2.97) 

(2.98) 

(2.99) 


4? = [(2* f + v n )U + (2 + ~ tr)] /a (2 >, 

Pzi d = [-3(1 - i/c)(l “ - (3 - u c - 2i/ ff )C„ + (1 - V()Kr ~ &)1 /* (2) > 

4? = [(*< + 2*/,)e, + (2 + *<)(*< - tr)] /*<*>, 


P^zz — [3(1 — ^c)(l + J'c + ^i?) + (3 + v c + 2*/,)^ + (1 — - £<)]/^ 2 >. (2.100) 


Then it can be shown that Eqs. (2.53)-(2.55) are equivalent to 


■S 1 * - 5 (4 ,, +tf ) +<?>). 

(2.101) 


(2.102) 

and 

(<ff /J =41 ) 4 ,) +pa ) 4 1) + 4 1 3 ) 4 1) , 

(2.103) 

where ( j,k,n + 1/2) € 0\. Also Eqs. (2.56)-(2.58) axe equivalent to 


= 1 (4 2) + 4 2) + 4 2) ) , 

(2.104) 

(*‘?« , = 4?4 J, +4?4 2) +4?4 2) . 

(2.105) 

and 

(“J )"fc 1 = P3 2) ^ /( l 2) + P3?4 2) + 4 2) 4 2) > 

(2.106) 


where (j, k, n + 1) € 0,2. Eqations (2.101)-(2.103), and Eqs. (2.104 )-(2.106), respectively, 
can be expressed as Eqs. (2.82) and (2.83) if one defines 


/ \ 


\ 


Q[ k) 






k as 1,2, and t = 1,2,3. 



(2.107) 


This follows from the fact that is the product of the row matrix which appears on 
the right side of Eq. (2.107) and one of the column matrices defined in Eq. (2.73). With 
the above definition, the marching procedure in the current a-fi scheme, i.e., the scheme 
which solves Eq. (2.1) with fi > 0, is formed by applying the marching steps defined by 
Eqs. (2.82) and (2.83) successively. Obviously, Eqs. (2.84) and (2.85) are also valid for the 
a-/i scheme. 
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Note that the a-y scheme can also be used to solve Eq. (2.1) with y = 0 as long as 
A^) ^ 0 and A^ / 0. According to Eqs. (2.32), (2.86) and (2.87), the last two conditions 
are reduced to Eq. (2.66) if y = 0. The current a scheme, however, is still applicable even 
if Eq. (2.66) is violated. 

The a-y scheme has several nontraditional feat tires. They are summarized in the 
following remarks: 

(a) Both the a and the a-y schemes have the simplest stencil in each of their two marching 
steps, i.e., a tetrahedron with a vertex at the upper time level and the other three 
vertices at the lower time level. 

(b) Each of the conservation conditions Eqs. (2.53)-(2.58) represents a relation among the 
marching variables associated with only two neighboring SEs. This is a fundamental 
difference between the current method and other traditional methods. 

(c) For both the a and the a-y schemes, we have 

q (j, k, n + 1) — ► q(j, k, n) as At — * 0, ( j , k, n ) £ ft, (2.108) 


if a, /x, and ax are held constant. The above property usually is not shared by other 
schemes which uses a mesh that is staggered in time, e.g., the Lax scheme [9, p.97]. 
The proof of Eq. (2.108) follows from Eqs. (2.84) and (2.85), and the fact that, as 
A f — ► 0, 


(d) 


(e) 


qPq? 


-* 0, 

ef'Qf - o, 

0, 

qWq® - 0, 

- * o? 

- o, 

(e^eS 2 ’ + + q^q?) 

->0, 

Q[ 2) Qi l> - o, 

ofoS 1 ’ - o, 

q< 2) q<” 0, 

— » 0, 

© 

T 

Sc 

o* 

/“V 

c* 

(oS 2 ) oS 1 ) +ofoi 1 ) +ofo| 1) ) 


—* lj 

-*• 1. 
(2.109) 


Alternatively, Eq. (2.108) can be proved using a line of argument involving flux con- 
servation similar to that given in the last paragraph on p.8 of [5]. 

Both the a and the a-y schemes are so called “two-way” marching schemes [5, p.ll]. 
In other words, the same flux conservation relations Eqs. (2.12) and (2.13) can be used 
to construct the backward time marching versions of the a and a-y schemes. More 
discussions on this subject are given in Appendix C. 

It will be shown in Sec. 5 that the a scheme is stable if 


|i^| < 1.5, |i/,| < 1.5, and \v$ + u v \ < 1.5. (2.110) 

As depicted in Fig. 9, the domain of stability defined by Eq. (2.110) is a hexagonal 
region in the v$-v v space. Moreover, it will be shown that (i) Eq. (2.110) can be 
interpreted as the requirement that the physical domain of dependence of Eq. (3.1) 
should fall within the numerical domain of dependence; and (ii) the a scheme is 
neutrally stable, i.e., free from numerical diffusion, if it is stable. 
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Also it will be shown that the a-p scheme is unconditionally stable for the pure 
diffusion case, i.e., when a = 0. The stability conditions of the a-p scheme with p ^ 0 
and a/0 are more complex. They also will be discussed in Sec. 5. 

Finally, it should be emphasized that, with the aid of Eqs. (2.19)-(2.22), (2.24), and 
(2.25), the a and the a-p schemes can also be expressed in terms of the marching variables 
and the coefficients tied to the coordinates ( x,y ). In other words, the coordinates 
are introduced solely for the purpose of simplifying the current development. The essence 
of the a and a-p schemes, and the schemes to be introduced in the following sections, is not 
dependent on the choice of the coordinates in terms of which these schemes are expressed. 
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3. The a-e Scheme 


The a scheme is neutrally stable and reversible in time. It is well known that a 
neutrally stable numerical analogue of 


du 

dt 


du vu. _ 
» r=0 ’ 


du 

dy 


(3.1) 


such as the a scheme, generally becomes unstable when it is extended to model the Euler 
equations. It is also obvious that a scheme that is reversible in time cannot model a physical 
problem that is irreversible in time, e.g., an inviscid flow problem involving shocks. In this 
section, we consider Eq. (3.1) and attempt to modify the a scheme such that it can be 
extended to model the Euler equations. 

To proceed, note that the CEs used in Sec. 2 will not be used in this section. As a 
result, Eqs. (2.12) and (2.13) will no longer be assumed. Instead, the CEs to be used are 


CE< 1) (j,k,n + 1/2) 

= [CES ,) (j,i,n + l/2)] U [C£$ 1) 0\*,n+l/2)] U [cE^O.i.n + 1/2)] , (3 ' 2) 


where (j, fc,n + 1/2) € fli, and 


C£ W (i,fc, " + !)=' [cEf’O'.i.n + l)] U [CJS^’O, *," + !)] U [CE< 2) y, Jb, n + 1)] , 

( 3 -3) 

where (j, k, n + 1) € We shall assume that the total flux leaving the boundary of any 
new CE vanishes, i.e., 


i 


S(CEW(j,k, n+l/2)) 


h*-ds = 0, 


and 


i 


h* -ds = 0. 


S(CE<s)(i,t, n +l)) 


(3.4) 


(3.5) 


Obviously, (i) E 3 can be filled with the new CEs, and (ii) the total flux leaving the boundary 
of any space-time region that is the union of any new CEs will also vanish. 

Moreover, it can be shown that Eqs. (3.4) and (3.5), respectively, are equivalent to 
Eqs. (2.67) and (2.70). 

Proof: By subtracting the expressions on the right sides of Eqs. (2.53)-(2.55), respectively, 
from those on the left sides, and then multiplying the results by 2wh[Z, we obtain the fluxes 
leaving CE^(j, k,n + 1/2), CE^O*, k,n + 1/2), and CE $P(j, k,n + 1/2), respectively (see 
Appendix A). Because the flux leaving an interface from the CE on one side is the negative 
of that leaving the same interface from the CE on the other side, it is easy to see that 
the flux leaving CE^(/, k,n- j- 1/2) (which is the union of CE^(j, k,n + 1/2), ^ = 1,2, 3) 
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is the sum of the fluxes leaving CE^(j, k,n + 1/2), £ = 1,2,3. Thus, the flux leaving 
CE (1) (j,Jfc , n + 1/2) can be obtained by subtracting the sum of the expressions on the right 
sides of Eqs. (2.53)-(2.55) from that on the left side, and then multiplying the result by 
2wh/3. Furthermore, we have 


(*)± (*)± (*)±_ 3 
°11 ' °21 + °31 — 


* = 1 , 2 , 


(3.6) 


and 


_(*)± , _(*)± (k)± _ (*)± (*)± (k)± _ , jo 

°12 T a 22 + °32 ~ °13 + "23 + a 33 ~ u , tC— 1,Z. 


(3.7) 


With the aid of the above considerations, and the fact that the expressions on the left sides 
of Eqs. (2.53)-(2.55) axe all evaluated at the same mesh point (j, k,n + 1/2), it becomes 
obvious that Eq. (3.4) is equivalent to the statement that is 1/3 of the sum on the 

right sides of Eqs. (2.53)-{2.55). By using Eqs. (2.35)-(2.43) and the assumption n = 0, 
we arrive at the conclusion that Eq. (3.4) is equivalent to Eq. (2.67). By invoking a similar 
argument involving Eqs. (2.56)-(2.58), and (2.44)-(2.52), we also conclude that Eq. (3.5) 
is equivalent to Eq. (2.70). QED. 

As a result, Eqs. (2.67) and (2.70) are shared by the a scheme and the current modified 
scheme. In this section we shall describe how the other equations in the a scheme i.e., 
Eqs. (2.68), (2.69), (2.71), and (2.72), can be modified such that the numerical diffusion 
of the resulting new scheme can be controlled by an adjustable parameter e. 

To proceed, for any (j, k,n + 1/2) € fii , let 


/ n+l/2 def f . At \ /«, n\ 

U j+l/3,*+l/3 ~ ( U +"5" U <) J (^-®) 

' \ Z / j+l/3,k+l/3 

i n+l/2 def f . At \ , . 

U j-2/3,k+l/3 = U + T Ut J ’ (3 - 9) 

\ * / i-2/3,*+l/3 

and 

(“ + f“0" ' (310) 

\ * J j+l/3,fc— 2/3 

By their definitions, u j+t /3 *+i/ 3 » u j- 2/3 *+i/ 3 ’ u j+t /3 *- 2/3 can considered as the 
finite-difference approximations of it at (j + l/3,fc + l/3, n + l/2), (j— 2/3, fc + 1/3, n + l/2), 
and (j + 1/3, k — 2/3, n-f 1/2), respectively. With the aid of Eqs. (2.10), (A. 10) and (2.31), 
Eqs. (3.8)-(3.10) imply that 


1 n+l/2 
u j+l/3,fc+l/3 





n 

j+1/3,*+1/3 ’ 


(3.11) 


u 


f n+l/2 
j-2/3,k+l/3 


it — 2 



+ u i u i) 


n 

i-2/3,fc+l/3 ’ 


(3.12) 
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and 


u 


t n+l/2 
j+l/3,fc-2/3 



+*v4) 


n 

j+l/3,k-2/3 ' 


(3.13) 


In Fig. 10(a), P, Q, and B are three points in the space. Using the coordinates 
given in the same figure, it can be shown that these points are on a plane represented by 


« = + K)£ 1/2 (-? - **,) + «'r t +i/2 , 


(3.14) 


where 

in+ 1/2 def 1 / fn+1/2 . ln+1/2 , fn+1/2 '\ 

M i,* — 3 ( U i+l/3,fc+l/3 + U j-2/3,k+l/3 + U >+l/3,fc-2/3j > 


(3.15) 


and 


K)”t ,/2 = 


( u > \ n + 1 / 2 4| f 

\ U T})j,k 


(' n+l/2 
\ U j+l/3,k+l/3 


i n+l/2 

^1+1/3, k+1/3 


«£&»/,) /AC. 


(3.16) 

(3.17) 


Equation (3.14) implies that point 0 depicted in Fig. 10(a) is also a point on the plane 
spanned by P, Q, and R (i.e., the plane that contains P, Q, R). Moreover, for every point 
on the plane represented by Eq. (3.14), including point O, we have 



(OS 1 *. 


and 



= {u 


I \ n+l/2 

v)j,k 


(3.18) 


As a result of the above considerations, and (uj,)?* ^ can be considered 

as the finite-difference approximations of tt, du/d(, and du/dr) at the mesh point (j, k,n + 
1/2), respectively. Note that u^ +1 ^ 2 generally is different from which is defined 

by Eq. (2.67). Furthermore, the former has no role in the future development. Let 


(.,•+ '\ n + 1 / 2 4*1 ^ ( u 1 \ n + 3/2 
l U C )j,k — Q ’ 


and 


( u '+\ n+1/2 d ± { ¥ L ( u ' r +1/2 

1 U 1 h,k Q \ U ri)j,k 1 


and 


(3.19) 


/,o+'» n + 1 / 2 1 / J 1 ) 

V u < 'i,* “ 3 V 1 



and 


f«° + V l+1 / 2 d — ^ (V 1 ) — 

is >J,k “ 3 V S 1 5 3 / • 


(3.20) 


where s^, s^, and are the expressions on the right sides of Eqs. (2.60)-(2.62), 
respectively. Because these expressions are functions of the marching variables at the nth 
time level, so are (u° + )J^ 1/ ^ 2 and (u®"*")”^ 1 ^ 2 . Moreover, Eqs. (2.68), (2.69) and (3.20) 

imply that the counterparts of (u^)"^ 1 ^ 2 and (u^)J^ 1 ^ 2 in the a scheme are (w^)”^ 1 ^ 2 
and (ujj + )"t 1/2 , respectively. 
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Let 


and 


(A,+);+ ,/2 =' 2 [(«' + )": I/2 - K + )"t 1/2 ] , (j,k,n + 1/2) £ fi„ 
«)£' /2 = 2 [« + )”t 1/2 - («; + )"t 1/2 ] . U,k,n + 1/2) € n,. 


(3.21) 

(3.22) 


Then, with the aid of Eqs. (3.11)-(3.13), (3.16), (3.17), (3.19), and (3.20), one has 


(d “< )>.‘ 1/2 = | [(“ + 4 “< ” 2U? )”- V 3, t+1 /3 “ (“ ~ 2U < “ 2 <)] 


j+ l/3,fc+l/3 


, (3.23) 


and 


(<K );,£ 1/2 = 5 [(«• - 2 “? + 4 <)" +1/3ii - 2/ 3 - - 2 “? - 2 <)J. 


7+l/3,*+l/3j 


. (3.24) 


Thus both (du + )”^ ly/2 and {du^)"^ 1 ^ 2 are functions of the marching variables at the nth 
time level. As will be shown shortly, they play a key role in the first marching step of the 
modified scheme. 

Next we consider Fig. 10(b). For any ( j , k,n + 1) G fit, let 


,„ +1 w a( \ n+,/2 

U >-l/3,Jfc-l/3 ~ l U “*" 0 Ut ) ’ 

\ * / ;-l/3, Jb-1/3 


(3.25) 


tt 


i n+l 
J+2/3,lr— 1/3 


def / , A f y 


and 


At y +1/2 

7+2/3,*-l/3 
At \ n+1/2 

,*+2/3 

By their definitions, u'”^ 3 Jt _ 1 ^ 3 , and u'”^ 3 j k+2 / 3 can be considered as the 

finite-difference approximations of u at ( j — 1/3, k - 1/3, n + 1), (j + 2/3, k — 1/3, n + 1), 
and (j — 1/3, k + 2/3, n + 1), respectively. With the aid of Eqs. (2.10), (A. 10) and (2.31), 
Eqs. (3.25)-(3.27) imply that 


d m( , At y +,/2 

“l-/3, i+2/ 3 - (“ + 2 “•) ,_ I/3 , 


(3.26) 


(3.27) 


u 


i n+l 
7+2/3, 


x-in+l/2 

+ v *) u t ) i 

^ ^ J >— 1/3,*— 1/3 

(3.28) 

. \ 1 "+1/2 

+ ^ W*-,/, ’ 

(3.29) 
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and 


“I-mVm-2/3 = [“ - 2 (‘W' + 

In Fig. 10(b), the points P, Q, and R are on a plane represented by 

“ = (“DS'K-W) + K)7.t’0» - i*>») + «J5 +1 . 

where 

/n+ l def 1./ /n+l ,,/n+l 4-M ,n+1 

_ 3 V U i-l/3,lfc-l/3 + U y+2/3,Jfc-l/3 + U j-l/3,k+2/3j > 


,fc+2/3 


(3.30) 


(3.31) 

(3.32) 


and 



/ n+l 

/+2/3,fc — 1/3 


/n+l 

j-l/3,k+2/3 


U 


/n+l 

j-l/3,k-l/3 



u' n+1 

U j-l/3,k-l/3 



Moreover, we introduce the current counterparts to Eqs. (3.19) and (3.20), i.e., 


(3.33) 

(3.34) 


(V + 'l n+1 — f — (V i n+1 
- 6 vVy,Jfc » 


“d (OS 1 = y«)S’. 


(3.35) 


and 


(« 


o+\n+l 

C A*,* 



and 


/'.,°+'\n+l def 



(3.36) 


where Sj 2 \ s^\ and are the expressions on the right sides of Eqs. (2.63)-(2.65), 
respectively. Because these expressions are functions of the marching variables at the 
( n + l/2)th time level, so are (tx£ + )”‘£ 1 and (w^ + )“ f * 1 - Moreover, Eqs. (2.71), (2.72) and 

(3.36) imply that the counterparts of (u^)"! 1 and (u^ + )”£ 1 in the a scheme are (u^)"^ 1 
and („;+)*« respectively. 

With the above preparations, the current counterparts to Eqs. (3.21) and (3.22) are 


and 



\n+l def 
)j,k ~ 



-(“C + )”,t 1 ]’ (>,*.»+ 1) 6 Sl 2 , 


(3.37) 


Mf)# 1 d 4 f 2 



“(“ 


o+\n+l 

1 '}* I ’ 


( j , k,n 1) € 


(3.38) 


respectively. With the aid of Eqs. (3.28)-(3.30) and (3.33)-(3.36), Eqs. (3.37) and (3.38) 
imply that 




[( 


tx + 2u 


+ 

C 


+ 



n+l/2 

j-l/3,k-l/3 


U - 4 U+ 



n+l/2 

i+2/3,Jfc — 1/3J ’ 


(3.39) 
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and 



n + 1 _ 1 
i,k 3 


[( 


u + 2u 


+ 

C 


+ 



n+l/2 

i-1/3,*— 1/3 


[u + 2u 


+ 

C 



n+l/2 

j- 1/3, *+2/3 


(3.40) 


Thus both (dtx^)”* 1 and (<fu+ are functions of the marching variables at the (n+l/2)th 
time level. 

The modified scheme can now be stated using the above definitions. It consists of two 
marching steps. The first is formed by Eq. (2.67), 


(OS* 


— fT/°+V ,+1 / 2 
)j,k 


+ e(du 


+ XT.+ 1/2 

< )j,k ’ 


(3.41) 


and 


(OS* = (OS* + «(<*OS* 


(3.42) 


where ( j,k,n + 1/2) E Cti , and e is an adjustable parameter. It was explained earlier that 
the expressions on the right sides of Eqs. (3.41) and (3.42) axe functions of the marching 
variables at the nth time level. Moreover, according to Eqs. (3.21) and (3.22), Eqs. (3.41) 
and (3.42) can also be expressed as 



(OS* = « + >S* + ( £ - l/2)(rf U +)“t 1/2 , 

(3.43) 

and 

<os* = « + )S* + (« - 1/2 wos* 

(3.44) 

respectively. 

The second marching step is formed by Eq. (2.70), 




(3.45) 

and 

(os* = (os 1 + ^os 1 - 

(3.46) 

where (j,k,n + 1) E ^ 2 - It was explained earlier that the expressions on the right sides 
of Eqs. (3.45) and (3.46) are functions of the marching variables at the (n + l/2)th time 
level. Furthermore, according to Eqs. (3.37) and (3.38), Eqs. (3.45) and (3.46) can also be 
expressed as 

(OS ’ + ( 347 ) 

and 

(OS 1 = K + )S* + ( £ - WOS 1 - 

(3.48) 


respectively. 
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Because the modified scheme is characterized by two parameters a and e, hereafter it 
is referred to as the a-e scheme. 

At this juncture, note that: 

(a) With the aid of Eqs. (3.20) and (3.36), it is seen that Eqs. (3.41), (3.42), (3.45), and 
(3.46), respectively, are reduced to Eqs. (2.68), (2.69), (2.71), and (2.72) when e = 0. 
As a result, the a-e scheme becomes the a scheme when e = 0. 

(b) For the special case with e = 1/2, Eqs. (3.43), (3.44), (3.47), and (3.48) are reduced 
to the forms that represent the finite-difference approximations defined in Eqs. (3.16), 
(3.17), (3.19), (3.33), (3.34) and (3.35). However, Eqs. (2.67) and (2.70), which are 
independent of e and therefore always part of the a-e scheme, are the results of the 
flux conservation conditions Eqs. (3.4) and (3.5). 

(c) With the aid of Eqs. (2.30) and (3.23), Eq. (3.41) can be rewritten as 


(“<) 


n+1/2 

i,* 



+ ) n+l/2 





n 

j-2/3,Jk+l/3 


f6u 2a v y 

\ A C A C / j+l/3,Jfc+l/3 

(3.49) 


Let (i) iiJ_ 2/3> * +1/ 3, K)"_ 2/3,k+i/3 “d («,)y_2/3,fc+i/3 be identified with the values 
of u, du/d( and du/drj at the mesh point (j — 2/3, k + 1/3, n), respectively; and (ii) 
u i+i/3,*+i/3> ( u c)”+i/3,fc+i/3 («»)"+i/3,fc+i/3 be identified with the values of u, 

du/d( and du/dr) at the mesh point (j + 1/3, A: + 1/3, n), respectively. Then it can 
be shown that the expression within the brackets on the right side of Eq. (3.49) is 
0(a£,at/). Furthermore, because Eq. (2.28) is applicable within SE (j,k,n) only, the 
expression that is enclosed within the first bracket on the right side of Eq. (2.28) is 
0(a(, At). From the above considerations, one concludes that the addition of the 
extra term involving e on the right side of Eq. (3.49) may result in errors that axe 
second order in a£, at/, and At. In other words, the addition of the term involving 
e does not result in a scheme of lower order of accuracy. A similar conclusion is also 
applicable to Eqs. (3.42), (3.45), and (3.46). 


(d) Equations (3.16), (3.18) and (3.19) imply that (u^)"^ 2 is proportional to the di- 
rectional derivative along the ^-direction on the plane spanned by points P, Q, and R 
depicted in Fig. 10(a). According to Eq. (3.21), (du*)"^ 1 ^ 2 is twice the difference be- 
tween (<4 + );t 1/2 and its counterpart in the a scheme. Note that the variable (du z )", 
that appears in Eqs. (3.2) and (3.10) of [5], plays a role in the 1-D a-e scheme [5] 
similar to that of (du^)"^ 1 ^ 2 in the present 2-D a-e scheme. It can be shown that 
(duz)'j is equal to the difference between two slopes. The first is the slope of a line 
spanned by 1/r2 and u'” ^ 2 , which are defined in Eq. (3.11) in [5]. The second is 
the counterpart of the first in the 1-D a scheme. Thus the 2-D a-e scheme is a natural 
extension of the 1-D a-e scheme. 
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The a-e scheme will take the form of Eqs. (2.82) and (2.83) if 


QW d 4 f - 
Vl 3 


fl-Vc-Vr, -(1 - - l/„)(l + U ( ) -(1 - I/ C - U v )( 1 + V n ) ^ 

1 — e -(1 + - 2e) -(1 + u n - 2e) 

V 1 — c -(l + ^-2e) — (1 + v v — 2e) / 


, (3.50) 


„<D a.. i 

V 2 3 


„(!> M 1 
3 


/ l + i/ c (1 + z/ c )(2 - i/ c ) -(l + ^)(l + ^)\ 
-(1 - e) -(2 - i/ c - 4e) 1 + v n - 2e 

\ 0 0 0 / 
/ H-I/ V -(1 + I/„)(H- I/ C ) (l + u v )(2-v v )\ 
0 0 0 
\ —(1 — e) 1 + i/£ — 2e —(2 — — 4e) / 


(3.51) 


(3.52) 


/ 1 + i/ c + i/. 


q ( 2) def 1 

3 


-(1 - e) 


(1 + i/ c + i/,)(l - i/ c ) (1 + i/ c + i/,)(1-i/,)\ 
-(l-i/ c -2e) -(1 - i/„ - 2e) 


V -(1 - e) -(1 - I/ C - 2e) -(1 -v v - 2e) / 


(3.53) 


and 


/l-v< 


n (2) def 1 

- 3 


1 — e 


-(1-^X2 + ^) (1 -z/ c )(l-i/,)\ 

-(2 + 1 /< - 4e) 1 - - 2e 


\ 0 0 


0 / 


/!-!/, (1 - i/,)(l - i/ c ) 


Qf d =' 5 


0 


0 


-(1 - i /„)(2 + I /,) \ 

0 


\ 1 — e 1 — u^ — 2e —(2 + i/, — 4e) / 


(3.54) 


(3.55) 


Note that: 

(a) The matrices defined above axe reduced to those defined in Eq. (2.80) when e = 0. 

(b) With the above definitions, Eqs. (2.84) and (2.85) are also valid for the a-e scheme. 
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In Sec. 5, it will be shown that (i) the a-e scheme is unstable if e < 0 or e > 1, and 
(ii) numerical diffusion increases as e increases, at least in the range of 0 < e < 0.65. In 
order to suppress numerical oscillations near a discontinuity, one may be forced to choose 
a large e. However, with such a choice, the smooth part of a solution may become highly 
diffusive. To solve this dilemma, in the following, we shall construct a generalization of 
the a-e scheme. 

To proceed, let (j, k,n + 1/2) 6 and consider Fig. 11(a). This figure is essentially 
identical to Fig. 10(a) except that point O in the latter is replaced by point O* in the 
former. The coordinates of point O* are (jaC, ^Ar/, u"^ 1 ^ 2 ) where is that defined 

in Eq. (2.67). Thus point O* generally is not on the plane spanned by points P, Q, and R. 
Let planes #1, #2, and #3, respectively, be the planes spanned by the following trios of 
points: (i) points O*, Q, and H; (ii) points O*, R, and P ; and (iii) points O*, P, and Q. 
Then in general these planes differ from each other and from the plane spanned by points 
P, Q , and R. In the following, first we shall study the former three planes. 


As a preliminary, let 


def m+ 1/2 n+1/2 

X 1 _ >fc+1 /3 w j,* > 

(3.56) 


def,„+l/2 _ n+1/2 

x 2 — U j_2/3,*+1/3 U j,k > 

(3.57) 


def l n+1/2 n+1/2 

*3 — u y+l/3,*r-2/3 U j,k » 

(3.58) 

= 

(«? ) )J» 1/, = -( 2 *»+*.)/aC, 

(3.59) 

i 

( u a ) ) n+ 1 / 2 M _ {XJ + 2x 3 )/&Tj, 

(3.60) 

= 

(«< ,) )S ,/1 = (2* .+*»)/aC, 

(3.61) 

- 

(«i 2) )7.t 1/2 = (*i - *»)/«(. 

(3.62) 

| 

(»f )m 1/2 = (*l - 

(3.63) 


and 



K 3) )M 1/2 = f ( 2*1 + x 2 )/Arj. 

(3.64) 


Moreover, let 


- 

(4'^r' 2 = (4 °)m 1/2 § + (“i' ) );,r /2 <= 1 - 2 , 3 , 

(3.65) 

- 


ia 

t 


30 


II i. IIUlllll lllkiil IBatUklili iii Uli LI, 



and 




dy 


£ = 1,2,3. 


(3.66) 


With the aid of Eqs. (2.20) and (2.22), we have 

r U ( / )l n+1/2 = ^ fti ( V +1/2 4 - ( u W ) n+1/2 

^ x *i’ k 2w ' U < )j ’ k + 2u> ^ 7 j ' k ’ 

and 


£ = 1,2,3, 


(3.67) 


/' t ,(*h"+ 1/ 2 _ (w + 6 )aC , (<Kn+l/2 (w-6)Aiy r fflvn+l/2 
ii.fc - 2wh + 2wh { " } j’ k ’ 

Combining Eqs. (3.59)-(3.64) with Eqs. (3.67) and (3.68), one has 


£ = 1 , 2 , 3 . 


(3.68) 


(4 1) )i 1 t 1/2 = -^(^2+X 3 ), 

(3.69) 

/ UK"+i/2 (36 + w)x 2 + (36 - u?)x 3 

> ** ~ 2 wh 

(3.70) 

(uW) n+1/2 - — 
>* k ~ 2w ’ 

(3.71) 

( f2Un+l/2 (36 + w)xi + 2WX 3 

{Uy } i’ k 2 wh ’ 

(3.72) 

fu< 3 h n+1/2 - 3X1 

(3.73) 

/. (3Un+l/2 ( w - 36)si + 2wx 2 
)j ’ k 2 wh 

(3.74) 


With the above preparations, it can be shown that planes #£, £ = 1,2,3, are repre- 
sented by 


» = (-W* (C - JAC) + (»<' > );,t J/2 (V - iA,) + u£ ,/2 , 6 = 1,2,3, (3.75) 

respectively, if the coordinates (£, 77) axe used. Alternatively, they can be represented by 


“ = (“i 0 )”^ 2 (* - *>,») + (V - «,.) + t = 1,2,3, (3.76) 
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respectively, if the coordinates ( x , y) are used. Using Eqs. (3.75) and (3.76), one concludes 
that, at any point on plane # £, £ = 1,2, 3, we have 


and 



(3.77) 


(3.78) 


Note that Eq. (3.77) is the current counterpart to Eq. (3.18) which is applicable to any 
point on the plane spanned by points P, Q, and R. Let Vu be the gradient of u. Then 
Eq. (3.78) implies that, at any point on plane # £, t = 1, 2, 3, we have 


|V»| = («,) 


n+1/2 def 
j,k 


+ ( 4 0 ) 2 


n+l/2 

y.fc 


(3.79) 


To proceed further, we introduce the current counterpart to Eq. (3.19), i.e., 

( u »+)"+> /2 ^ (tl ('))"+'/ 2 , and (u<« + );+ 1/2 = ^K 0 )“,r /2 . 

Then Eqs. (3.16), (3.17), (3.19), and (3.56)-(3.64) imply that 

(«?)S w *-i[«r+-r+-r]^ 

and 

(“', + )"t i/2 = l [4 i>+ + «? >+ + «? )+ ] 2 

i.e., (i) (u' + )”£ 1/2 is the simple average of 

/ (i)+xn+l/2 / (2)+xn+l/2 

'j,* > 

and (ii) (u'+ )"”£ 1/f2 is the simple average of 


| n+l/2 

,* 

n+l/2 

,k 


and 


,(3)+^n+l/2 
\k » 




and (n< 3 > + )”,: ,/2 . 


(3.80) 


(3.81) 


(3.82) 


(3.83) 


(3.84) 


The first marching step of the generalized a-e scheme will be formed using Eqs. (2.67), 
(3.43), and (3.44) except that (i) (u^ + )”^ 1 ^ 2 in Eq. (3.43) is replaced by the weighted 

average of those given in Eq. (3.83); and (ii) (uj* in Eq. (3.44) is replaced by the 
weighted average of those given in Eq. (3.84). The design of these weighted averages will 
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be guided by the requirement that the weight assigned to a quantity associated with plane 
is bigger if is smaller. This requirement is similar to that put forward on 

p.28 of [5]. 

For any a > 0, the weighted-average counterparts to (u ^)”* 1 ^ 2 and are 


f0, 


u 


nj-f- def 


= < 


and 


uj-f def 
“t = 


( + ( W«ff )+ + ( 

( #1#2 )° + (#2#3)“ + (#3#l)“ 


0 , 

(W*4 1H +(w4 2H +(t>it>2)°«i 3,+ 

(«1« 2 )“ +(«2«3>” + (Wl)“ 


if &i — 82 = #3 — 0; 
otherwise, 

if 9 i = $2 = $3 = 0 ; 
otherwise, 


(3.85) 


(3.86) 


respectively. Here, for simplicity, we suppress the indices j, k, and n + 1/2 which are 
associated with all symbols in Eqs. (3.85) and (3.86). Because the denominators of the 
fractions on the right sides of Eqs. (3.85) and (3.86) vanish if a > 0 , and any two of # 1 , 62 , 
and 63 vanish, consistency of the above definitions requires the proof of the proposition: 
#1 = #2 = #3 = 0 , if any two of $i, 9 2 , and $3 vanish. 

Proof: As an example, let 9\ = 9 2 = 0. Then Eq. (3.79) implies that = 0, 

t — 1,2. In turn, Eqs. (3.69)-(3.72) imply that xi = x 2 = X 3 = 0. # 3=0 now follows 
from Eqs. (3.73), (3.74), and (3.79). QED. 

Next we consider several special cases of Eq. (3.85). We have 


</ c i)+ , if », = 0 , 

0 

A 

Ci 

and # 3 > 0 ; 

u' 2)+ , if « 2 = 0 , 

O 

A 

H 

and #3 > 0 ; 

u f + , if = 0 , 

#! >0, 

and #2 > 0 . 


(3.87) 


Assuming #/ > 0, £ = 1,2,3, we have 


„„ (i/o,)^ 1 * + (VhYuf* + (i/o 3 )°u 

C (l/«i)“ +(l/ 0 2 )“ +(l/ 0 3 )“ 


(3.88) 


Thus the weight assigned to u^ + is proportional to (1 /#/)". By using Eqs. (3.81), (3.85), 
and (3.88), one arrives at the conclusion that 

u? + = u' + , if #! = # 2 = # 3 . (3.89) 
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Obviously Eqs. (3.87)-(3.89) Eire still valid if each symbol C in them is replaced by the 
symbol r/. 

On the smooth part of a solution, 6 \ , $ 2 , and $z axe nearly equal. Thus the weighted 
averages u“ + and u“ + are nearly equal to the simple averages u^ + , and u' n + , respectively 
(see Eqs. (3.81) and (3.82)). In other words, the effect of weighted-averaging generally is 
not discernible on the smooth part of a solution. 

Next let (j, k,n + 1) £ and consider Fig. 11(b). The third coordinate u”^ 1 of point 
O * is that defined in Eq. (2.70). Let planes #1, #2, and #3, respectively, be the planes 
spanned by the following trios of points: (i) points O*, Q, and R; (ii) points O*, R, and 
P; and (iii) points O*, P, and Q. In the following, we shall study these planes. 


As a preliminary, let 



-« n+1 

V 1 — U j-l/3,k-l/3 U j,k ’ 


(3.90) 


y 2 *“ U j+2/3,k-l/3 U j,k ■> 


(3.91) 


V3 = u'”|/ 3) fc + 2/3 “ 


(3.92) 


(“C”® 1 =(2»2 + »3)/AC, 


(3.93) 


=* ( V2 +2y 3 )/Aij, 


(3.94) 


(“C 2 ’)?! 1 = -(2».+»3)/a<, 


(3.95) 


(“i 2) )”,t I =* (vs - yijM, 


(3.96) 




(3.97) 

and 



(3.98) 

Moreover, let 

= K' ) )"i , § +«'>)?, r|> 

£=1,2,3, 

(3.99) 

and 


£ = 1,2,3. 

(3.100) 
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With the above definitions, Eqs. (3.67) and (3.68) remain valid if each upper index n + 1/2 
is replaced by n + 1. As a result, Eqs. (3.93)-(3.98) imply that 


(■w -£<»+*). 

(3.101) 

/„.(i)\n+i (36 + w)y 2 + (36 - w)y 3 

>>,t ~ 2wk 

(3.102) 

fu< 2 )r +1 - 3yi 

(3.103) 

/J2)\n+i (36 + w) yi + 2u>y 3 

“ 2 wh 

(3.104) 

( ( 3 Un+l _ 3 1/1 

( * )j ’ k ~ 2w ’ 

(3.105) 

/..rs^n+1 (w - Zb)yx + 2wy 2 

[U y >i’ k ~ ~ 2wh 

(3.106) 


With the above preparations, the earlier developments that involve Eqs. (3.75)-(3.89) 
can be repeated for the current case with the only change being the replacement of each 
upper index n + 1/2 by n + 1. Particularly, (u^ + )J^ 1 and («JJ' + )" ( * 1 can be defined using 
Eqs. (3.85) and (3.86) with the understanding that each symbol in these equations is 
associated with the mesh point (j,k,n + 1). 

The generalized a-e scheme, referred to as the weighted-average a-e scheme, can now 
be stated. It consists of two marching steps. The first is formed by Eq. (2.67), 



(“c )m‘ /2 = («r>M /2 + < e - l/2)(d»J)”t 1/2 . 

(3.107) 

and 

«)m 1/2 - wr)5 ,/2 + < e - 1/2)(*4 )m I/2 > 

(3.108) 

where (/, k, n 

+ 1/2) E fii. The second is formed by Eq. (2.70), 



(“< )"t‘ = («r f )5 1 +(«- 1/2KJU+)”; 1 , 

(3.109) 

and 

tost 1 = (-r )'T +(<- i/zxoit 1 . 

(3.110) 

where (/, k, n 

+ 1) € Si 2 ’ 


Note that, according to Eq. (3.79), the evaluation of (9i) a does not involve a fractional 
power if a is an even integer. Because a fractional power is costly to evaluate, the use of 
the generalized a-e scheme is less costly when a is an even integer. 
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4. The Euler Solver 


We consider a dimensionless form of the 2-D unsteady Euler equations of a perfect 
gas. Let p, u, v, p, and 7 be the mass density, ^-velocity component, y - velocity component, 
static pressure, and constant specific heat ratio, respectively. Let 

u 1 = p, u 2 = pu, u 3 = pu, = p/( 7 - 1) + PO 2 + u 2 )/2, (4.1) 


SI = « 2 , 

SI = (7 - 1)«4 + (3 - 7 )(u 2 ) 2 /( 2 ui) - (7 - 1)(u 3 ) 2 /(2u!), 

SI [ = «2«3/«l, 

S 4 = 7«2U4/«1 - (l/2)(7 - 1)«2 [(«2) 2 + («s) 2 ] /(«l) 2 , 

/? = «3, 

S 2 = «2U3/«1, 

S 3 = (7 - 1)«4 + (3 - 7 )(u 3 ) 2 /( 2 u!) - (7 - 1)(u 2 ) 2 /(2u 1 ), 

and 

/? = 7 U 3 u 4 /ui - (l/2)(7 - l)u 3 [(u 2 ) 2 + (u 3 ) 2 ] /(ui) 2 . 
Then the Euler equations can be expressed as 


= 

dt dx dy 


m = 1, 2, 3,4. 


(4.2) 

(4.3) 

(4.4) 

(4.5) 

(4.6) 

(4.7) 

(4.8) 

(4.9) 


(4.10) 


The integral form of Eq. (4.10) in space-time E 3 is 

<p h m • ds = 0, m = l,2, 3, 4, (4.11) 

JS(V) 

where 

£.-(£./*,«»). m = 1,2, 3, 4, (4.12) 

are the space-time mass, x-momemtum component, y-momemtum component, and energy 
current density vectors, respectively. 

As a preliminary, let 


=' drjdut, and = dfl/du,, m,/ = 1,2, 3,4, (4.13) 


- def , 

■Ut = u t /u 1 , 


£ = 2,3,4, 


(4.14) 
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and 


u 2 d = (u 2 ) 2 + (u 3 ) 2 . (4.15) 

Let F x and F y denote the matrices formed by /£ t and t , m,£ = 1,2, 3, 4, respectively. 
Then 


/ 0 

- u 2 -( u 2 ) 2 


F x = 


0 


0 \ 


(3 - t)w 2 


-« 2«3 


«3 


(1-7)u 3 7-1 

u 2 0 


(7 - l)u 2 u 2 - 7« 2 u 4 7«4 - ^ 2 ■■ ■ [2 (u 2 ) 2 + « 2 ] (1 - 7)6263 762 J 


and 


(4.16) 


/ 


F y = 


-U 2 Uz 


7 — 1 -2 \ 2 

— s— “ - («ar 


0 

«3 

(1 - 7)«2 


1 

u 2 


(3 - 7)63 

7-1 


0 \ 
0 

7-1 


\( y - 1 )u 3 « 2 — 7W3M4 (1 — 7 )u 2 « 3 764 - — [2(« 3 ) 2 + u 2 ] 763 J 

(4.17) 

Because /£ and /^, m = 1,2, 3,4, are homogeneous functions of degree 1 [23] in «i, 
U2, «3, and u 4 , we have 


Sm = F, “<> “d fm=Yl 


(4.18) 


*=1 


/=1 


For any (x,y,f) G SE(j,fc,n), u m (x,y,t), /*(x,y,t), /»(x,y,t), and £ m (x,y,*), re- 
spectively, are approximated by «;(x,y,< ;j,*,n), /**(x,y,i ; j, fc,n), /^,*(x,y,t fc,n), 
— ♦ 

and h* m (x,y,t\j,k,n). They will be defined shortly. Let 

«m(*.y»*;i>M) d = (««)",* + («»»*)",*(* - ***) + («my)",*(y - y;,*) ^ 

+ («mt)",t(< - < n ), m = 1, 2, 3,4, 

where (u m )j k , (u mx )j k , (u mj ,)” >jt , and (u m <)" ( * are constants in SE(j,A;,n). Obviously, 
they can be considered as the numerical analogues of the values of u m , du m /dx, du m /dy , 
and du m /dt at (xj,yfe,t n ), respectively. 
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i ii III in 


m = 1,2, 3, 4. Let 


and 


Because (i) 



denote the values of /^, /^, 

fm,V 811(1 

r , when u m , m = 1,2, 3, 4, respectively, assume the values of (u m )" fc , 

(/».)”* = izuuhMh, 

t = i 

m = 1,2, 3, 4, 

(4.20) 

= 'EUZJhMl* 
1=1 

m = 1,2, 3,4, 

(4.21) 

= y tuZ,,)lkMlk, 
1=1 

m = 1,2, 3, 4, 

(4.22) 

(ft.)?.* = 

i= 1 

m = 1, 2, 3, 4, 

(4.23) 

c/i,);.. = 

i=i 

m = 1,2, 3, 4, 

(4.24) 

/=i 

m = 1,2, 3,4. 

(4.25) 

5 /m V- f x 

a» Jm ' e dx ’ 

m = 1,2, 3, 4; 

(4.26) 


and (ii) the expression on the right side of Eq. (4.20) is the numerical analogue of that on the 
right side of Eq. (4.26) at ( Xj,y*,t n ), can considered as the numerical analogue 

of the value of dfa/dx at (iy,y fc ,t n ). Similarly, (/£,)$*, (/£<)",*> (/&*)",*> (/&*)",*> 311(1 
(/mt)> jt can be considered as the numerical analogues of the values of df^/dy, df^/dt, 
dfj^/ox, df^/dy, and dfj^/dt at (xj,yk,t n ), respectively. As a result, we assume that 


t,n) = (/’)",* + (/;,);*(* - */,») + (/A, )",.(» - «,») 

+ (/i.)S »(<-<•), m = l,2,3,4, 


(4.27) 


and 


n) « (/»)»„ + (/',)",,(* - + (/»,)* t (y - y M ) 
+ (/i, )?,*(»- <”), m = l,2,3,4, 


(4.28) 
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Also, as an analogue to Eq. (4.12), we assume that 

h* m (x, y, t]j, k, n) = f (/”(*> y,t;j,k , n), f%*(x, y, t ; j, k , n), 


«m(*» »»*;;'»*»»))» 


m = 1,2, 3, 4. 


(4.29) 


Note that, by their definitions: (i) /^, /^, and t , m,£ = 1,2, 3,4, are functions 

(«m)j,*> ™ = 1,2, 3, 4; (ii) (/* x )£* and (/» *’)£*, m = 1,2,3, 4, are functions of ( u m )£* 
and (« mi )” t , m = 1,2, 3, 4; (iii) (/*,,)£* and (/&,,)£*, m = 1,2, 3, 4, are functions of 
(«*.)£* and (u my )l k , m = 1,2, 3, 4; and (iv) (/*,)£* and (/" ,)£* are functions of (u m )£* 
and (umt)” (fc , m = 1,2, 3, 4. 

Moreover, we assume that, for any (x, y, f) € SE(j, k, n), and m = 1, 2, 3, 4, 


d u *m (*> y J * ; h M) + df^(x,y,t;j,k,n) + df£(x,y,t;j,k,n) _ Q 


at 


ax 


ay 


(4.30) 


Note that Eq. (4.30) is the numerical analogue of Eq. (4.10). With the aid of Eqs. (4.19), 
(4.27), (4.28), (4.20), and (4.24), Eq. (4.30) implies that, for m = 1,2, 3, 4, 

= -uz,)h - — £ f r mf uu + « J " . (4.3i) 

<=i J ’ 

Thus (u m( )" ilt axe functions of («m)",t, («rai)"t, and (u my )" >Jk . From this result and the 
facts stated following Eq. (4.29), one concludes that the only independent discrete variables 
needed to be solved in the current marching scheme are (u m )" fc , (u mx )] k , and (u my )’- k . 

Consider the conservation elements depicted in Figs. 5(a) and 6(a). The Euler coun- 
terparts to Eqs. (2.12) and (2.13), respectively, are (i) 


/ 


' S(CE^(j,k,n+l/2)) 
where (j, k,n + 1/2) € ft^ and (ii) 


h* m -ds = 0, £ = 1,2,3, m = 1,2, 3,4, 


(4.32) 


/ 

Jst 


S(CE?\j,k,n+ 1 )) 


h* m -ds = 0, £ = 1,2,3, m = 1,2,3, 4, 


(4.33) 


where (j, k, n + 1) € ft 2 . 

Next we shall introduce the Euler counterparts of Eqs. (2.24), (2.25), (2.30), and 
(2.31). For any (j, k, n ) 6 ft, let 



u^)i t 


m,£ = 1, 2,3,4, 


(4.34) 
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and 



m,£ = 1, 2,3,4, 


(4.35) 


The normalized counterparts of those parameters defined in Eqs. (4.34) and (4.35) axe 

(/&)?,* = “d (/’+)”.* = (4.36) 


and 


Md ( 4 - 37 ^ 


In the following development, for simplicity, we may strip from every variable in an 
equation its indices j, k, and n if all variables are associated with the same mesh point 
(j, fc, n) G 0. Let F < > + and F n+ , respectively, denote the matrices formed by and 
— 1,2,3, 4. Let I be the 4x4 identity matrix. Then the current counterparts 
to Eqs. (2.35)-(2.52) are 


s (l } ± def J _ F C+ _ F n+, 

(4.38) 

£g ):t d = ±(I - F c+ - F'+)(I + F<+), 

(4.39) 

s (l)± def ± (j _ F <+ _ F n+)(7 + F” + ), 

(4.40) 

s£ )±d = I + F<+, 

(4.41) 

E^ )d: d = t(I+ F c+ )(27 - F c+ ), 

(4.42) 

s (i)± def + F C+)( i + F n+) ? 

(4.43) 

E^ ):t d = 7 + F r,+ , 

(4.44) 

sg ):t d = ±(7 + F*+)(7 + F<+), 

(4.45) 

E^ ):t d = t( 7 + F’ f+ )(2I - F T?+ ), 

(4.46) 

d = 7 + F<+ + F’ 7+ , 

(4.47) 

Eg ):t = f =F (7 + F c+ + F*+)(7 - F< + ), 

(4.48) 

Eg )d: d = q=( J + F< + + F V+ )(I - F v+ ), 

(4.49) 

s (2)± drf j_ F <+ 5 

(4.50) 

sg }± d = ±(I - F c+ )(27 + F c+ ), 

(4.51) 
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Eg* =‘=F(/-F< + )(/ -?>+), 

(4.52) 

e‘ 2 i ,± =' I -F1+, 

(4.53) 

Eg >±d = T(7-f’ + )(/-f< + ), 

(4.54) 

eg 1 * =' ±(/ - F”^)(2I + F»+). 

(4.55) 


and 


Note that, for the case /z = 0, Eqs. (2.35)-(2.52) become Eqs. (4.38)-(4.55), respectively, 
if (i) 1, i/{, and u v , are replaced by J, and F n+ , respectively, and (ii) and 

are replaced by E^* and E^^, m,£= 1,2,3, respectively. 

Equations (4.32) and (4.33) are evaluated in Appendix B. Let u, , and u + , respec- 

tively, be the 4x1 column matrices formed by u m , tz*£, and m — 1,2, 3, 4. Then the 
results can be expressed as: 

(a) Eq. (4.32) with i — 1: 


[s ( ,‘. ,+ ff + E ( i 2 )+ \ + + 2< iV + 5, + ] " t = [eS’/'*? + eS'/’s/ + +] " 

(b) Eq. (4.32) with i = 2: 

e£ h s+ Eg )+ ff+ + s‘‘ )+ tr+] “J 1/a = [s<V-ff+ sg’-ff/ + 


j+1/3,*+1/3 

(4.56) 


(c) Eq. (4.32) with i = 3: 


H-i/3 

(4.57) 


e£ )+ < 7 + Eg )+ <r+ + E< 1 3 )+ ff+]^ 1/2 = [s<;>-ir+ + E<‘>-ff+]” 


(d) Eq. (4.33) with £=\: 

■>( 2 )+.v+ _i_ v( 2 )+.-r+l n+1 


s ii ,+ff + sS? + ff < + + sS 2 3 )+ ff, + l" = [sf 1 , "tr+ Eg )- <r+ + Eg ) -ff+' 


(e) Eq. (4.33) with l = 2: 

[Eg )+ 5+Eg )+ u+ + E< 2 3 )+ iT+]" +1 = [E< 2) -ff+ S< 2) -ff+ + Eg )_ u+]” +1/2 

(f) Eq. (4.33) with t = 3: 


i+l/3,*-2/3 

(4.58) 


n+1/2 

j-1/3,*— 1/3 ' 

(4.59) 


j+2/3,*-l/3 

(4.60) 




v(2)+£ . E (2 )+ u + 4- £^ + ?7 + 

^31 u T 2j 32 I ^33 u i} 


1 n+1 

}>k 


S31 “ + E32 "iz^ -f E 
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n+ 1/2 

j-l/3,fe+2/3 ’ 

(4.61) 
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Here ( j , k, n+1/2) € Qi is assumed in Eqs. (4.56)-(4.58) while (j, fc,n + l)€ Oj is assumed 
in Eqs. (4.59)-(4.61). Note that the forms of Eqs. (4.56)-(4.61) are very similar to those 
of Eqs. (2.53)-(2.58). 

As a result of Eqs. (4.38)-(4.55), we have 


Eff* + s£ ):k + S^ ):fc = SI, k = 1, 2, (4.62) 

and 

si?* + Eg* + Eg’* = sfi* + Eg’* + Eg'* = 0, 1=1,2. (4.63) 

Equations (4.62) and (4.63) are the Euler counterparts of Eqs. (3.6) and (3.7), respectively. 
By summing over Eqs. (4.56)-(4.58), and using Eqs. (4.62) and (4.63) with k = 1, one 
concludes that, for any (j, fc, n + 1/2) G f2j, 


“ ir /2 


i{ E<+«t + s S , 3 , 'N + ]; +1/ 3, t+1/3 

+ [ e £>-* + e ®-«? + 4 V -<]*. 2/3il+1/3 
+ [eS>-* ■ + e£>-<? + 4 V - N + ]’ +1/3ii . 2/3 }. 


(4.64) 


As a result, can be evaluated in terms of the marching variables at the nth time 

level. Similarly, by summing over Eqs. (4.59)-(4.61), and using Eqs. (4.62) and (4.63) with 
k = 2, one concludes that, for any (j,k,n + 1) € fi 2 j 


“It 1 




+ S ( 1 2) "u+ + E 


( 2 ) jt+ 
13 “ij 


n+1/2 

j-l/3,k-l/3 


+ + Eg )_ u+ + Eg ) -tx+] " +1/2 


i+2/3,fc— 1/3 


(4.65) 


( s » 


As a result, can be evaluated in terms of the marching variables at the (n + l/2)th 
time level. 

For any (j, k,n + 1/2) € fti, the matrices (E^ + )J j ^ 1 ^ 2 , m,£ = 1,2,3, are functions 
of u"* 1/2 . Thus they are also functions of the marching variables at the nth time level. 

Assuming the existence of the inverse of each of the matrices (£^{ + )"^ x ^ 2 , m = 1,2,3, 
one can define 


gW 


-imrn 


a . yt 1 ) - f?+l ” 

l> u u + 2j 12 u* +- 2j 13 u v 


J/+l/3,*+l/3 ’ 


(4.66) 
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(4.67) 


5 2 0) = [(Sg )+ )”J /2 ] ' [s<i ) “i? + S^ ) 'ff < + + S< , , ) -5+]”_ 2/3> 


k+ 1/3 


and 


s, (,) = [(si; >+ )" 2 1/2 ] 1 [s< 1 . ) *ff + $$-*? + 


i+l/3,*-2/3 


(4.68) 


where the inverse of a matrix A is denoted by [A] 1 . As a result of their definitions, S^\ 
t = 1,2,3, can be evaluated using the marching variables at the nth time level. 

For any ( j,k,n + 1) G fyj, the matrices (E^ + )"* 1 , m, i — 1,2,3, are functions of 
• Thus they are also functions of the marching variables at the (n + l/2)th time level. 
Assuming the existence of the inverse of each of the matrices (Sm/)"! 1 , m = 1, 2,3, one 


can define 


5™ d =' [(s ( ,? + )” 


n+1 

k 


gW d _£ f 


-1 r 


= (4i H ) "l' [sg )- * + + s 23 )_a 


,(2) "u + sS?'«+ + S 


(2)- 

13 


^+l n + 1 /2 
n Jj-X/3, 


*-1/3 


-- + l n+1 / 2 


J j+2/3,*-l/3 


(4.69) 

(4.70) 


and 


3> W = [( 2 3? + )”^] 1 [4?~ff+sg ) '\ + + 4l > ~S + 


_ _j_l n +X/2 

j-l/3,*+2/3 


(4.71) 


As a result of their definitions, \ i = 1,2,3, can be evaluated using the marching 
variables at the (n + l/2)th time level. 

Using Eqs. (4.38), (4.41), (4.44), (4.47), (4.50), (4.53), and (4.66)-(4.71), Eqs. (4.64) 
and (4.65) can be recast as 

5“t 1/2 = 5 [( J - ■ F<+ - ■P”’ + )”t 1/2 5, (1) + (/ + F<+)"+ 1/2 5 2 (1) + (/ + F’ + )"j 1/2 5<‘>] , 

(4.72) 

and 


+ F <+ + 3™ + (I- •F’ <+ )”t 1 S} 2) + (I - F^ 1 3 


C(2) 

3 


(4- 73 ) 

respectively. Equations (4.72) and (4.73) are the Euler counterparts of Eqs. (2.67) and 
(2.70), respectively. 

Furthermore, Eqs. (4.38)-(4.55) imply that: 
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(a) At any (j, k,n + 1/2) € fti, we have 



and 

(b) At any (j, k, n + 1) € 


and 


Eg*] 1 Eg* = ± (7 + F<+) , 

(4.74) 

Eg*] ' Eg* = ± (7 + F'i+) , 

(4.75) 

E<'*] Eg* = T (2/ - F< + ) , 

(4.76) 

[E<V + ‘]' 1 Eg* = ±(/ + F’ + ), 

(4.77) 

Eg*] Eg* = ± (/ + F <+ ) , 

(4.78) 

[e' 1 1 ,+ ]”‘e' 1 * = t(2/-F’ + ). 

(4.79) 

2 , we have 


[sg*] Eg* = T (/ - 7*+) , 

(4.80) 

[e< 2 *] _1 e < 1 2 * = t (J'-f ,+ ), 

(4.81) 

[Eg*pEg*=±(27 + F< + ), 

(4.82) 

4i*]’ i 4 2 * = t(/-f , ' + ), 

(4.83) 

[Eg*] Eg* = T (/ - f*) , 

(4.84) 

[sg*] Eg* = ± (2 7 + F*) . 

(4.85) 
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By using Eqs. (4.74)-(4.85), Eqs. (4.56)-(4.61) imply that 


[«+(/ + F<+) u + + (7 + F v+ ) «+] ”^ /2 = 5 X (1) , 

(4.86) 

[«- (27 — 7 1< * + ) iT+ + (7 + F n+ ) iT+]“* 1/2 = S 2 (1) , 

(4.87) 

[« + (7+ F<+) u + - (27 - F" + ) tT+1 ” +1/2 = S 3 (1) , 

L S J 3,k 

(4.88) 

[s _ (/ _ j*+) s+ -(I- F'+) ff+1 = £<*>, 

(4.89) 

[u + (21 + F< + ) u+ - (I - F" + ) ff+1 " +1 = S 2 (2) , 

L J j y k 

(4.90) 

U- (I- F< + ) S+ + (21 + F’+) a+1 = 5 3 (2) , 

^ J 

(4.91) 


where ( j,k,n + 1/2) 6 fii is assumed in Eqs. (4.86)-(4.88), while ( j,k,n + 1 ) 6 Q 2 is 
assumed in Eqs. (4.89)-(4.91). Because s^\ and denote the 

expressions on the right sides of Eqs. (2.60)-(2.65), respectively, a comparison between 
these equations and Eqs. (4.86)-(4.91) reveals that the latter are the Euler counterparts 
of the former, respectively. 

Note that, by multiplying Eqs. (4.86)-(4.88) from left with (7 — F^ + — F T,+ )^ 1 ^ 2 , 

(7+F ( * + )”^' 1 ^ 2 , and (I+F v+ )^ 1 ^ 2 , respectively, and summing over the resulting equations, 
one can obtain Eq. (4.72). Similarly, one can obtain Eq. (4.73) from Eqs. (4.89)-(4.91). 
Moreover, (i) Eqs. (4.86)-(4.88) also imply that 


= 5 - sT) , (4.92) 

where ( j,k,n + 1/2) € fir, and (ii) Eqs. (4.89)-(4.91) also imply that 

(s< + )”/ = l$ 2) - S t (2> ) , rnd (u+)”*‘ = I ($> <2) - S™) , (4.93) 


where (j, k,n + 1 ) 6 fi 2 . By using the above results, and directly substituting Eqs. (4.72), 
(4.73), (4.92), and (4.93) into Eqs. (4.86)-(4.91), one concludes that: (i) Eqs. (4.72) and 
(4.92) are equivalent to Eqs. (4.86)-(4.88); and (ii) Eqs. (4.73) and (4.93) are equivalent 
to Eqs. (4.89)-(4.91). 

With the above preparations, an Euler solver can now be defined. It consists of two 
marching steps. The first is formed by Eqs. (4.64) and (4.92), while the second is formed 
by Eqs. (4.65) and (4.93). As explained earlier, S} k \ k = 1,2, and t = 1,2,3, become 
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known after and are evaluated using Eqs. (4.64) and (4.65), respectively. This 

Euler solver has a two-way marching nature similar to that of the a scheme. As a result, it 
must be neutrally stable, (i.e., no numerical diffusion) if it is stable. Because it is reversible 
in time, this solver cannot model a physical problem that is irreversible in time, e.g., an 
inviscid flow problem involving shocks. Hereafter, this new Euler solver will be referred to 
as the Euler a scheme. 

At this juncture, note that the Euler a scheme is greatly simplified by the fact that 
and u"! 1 , respectively, can be directly evaluated in terms of the marching variables 
at the nth and (n + l/2)th time levels (see Eqs. (4.64) and (4.65)). As a result, the 
matrices (E^ + )J^ 1 ^ 2 and (E^| + )”j[‘ 1 , which are nonlinear functions of and u?^ 1 , 

respectively, can be evaluated easily. In other words, nonlinearity of the above matrix 
functions does not cause a particular problem for the Euler a scheme. 

To explain how Eqs. (4.64) and (4.65) arise, note that 


£ 


S(CE(i>(j\*,n+l/2)) 


h* m • & - (h Tl + 1/2) 6 fti, 


(4.94) 


and 


£ 


S(CEW{j,k,n+l)) 


h*m‘ds = 0 . 


(j',fc,n + 1) G fi 2 , 


(4.95) 


respectively, are the direct results of Eqs. (4.32) and (4.33), the basic assumptions of the 
Euler a scheme. According to Eq. (3.2), CE^(j, k,n + 1/2) is the hexagonal cylinder 
A'B'C'D'E'F 1 AB CDEF depicted in Fig. 5(a). Except for the top face A'B'C'D'E'F ' , 
the other boundaries of this cylinder are the subsets of three solution elements at the 
nth time level. Thus, for any m = 1,2, 3, 4, the flux of h leaving CE^\j, k, n + 1/2) 
through all the boundaries except the top face can be evaluated in terms of the marching 
variables at the nth time level. On the other hand, because the top face is a subset of 
SE(j, fc, n - f 1/2), the flux leaving there is a function of the marching variables associated 
with SE (j,k,n + 1/2). Furthermore, because the outward normal to the top face has no 
spatial component, Eq. (4.29) implies that the total flux of h leaving CE^(j, k,n + 1/2) 
through the top face is the surface integration of u£, over the top face. Because the center 
of SE(j,k,n + 1/2) coincides with the center of the top face , it is easy to see that the 
first-order terms in Eqs. (4.19) do not contribute to the total flux leaving the top face. It 
follows that the total flux leaving the top face is a function of (u m )"^ 2 only. As a result 
of the above considerations, can be determined in terms of the marching variables 

at the nth time level by using Eq. (4.94) only. Similarly, if*** 1 can be determined in terms 
of the marching variables at the (n + l/2)th time level by using Eq. (4.95) only. Eqs. (4.64) 
and (4.65) are the direct results of Eqs. (4.94) and (4.95), respectively. 

In an extension currently under development, the mesh used is not uniform in space. 
As a result, point G' depicted in Fig. 5(a) generally is not the center of the top face referred 
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to earlier. To simplify the development, we have moved the center of SE(j, k,n + 1/2) to 
the center of the top face, i.e., away from point G'. 


Next we shall construct the Euler a-e scheme, i.e., the Euler version of the a-e scheme. 
For this scheme, we shall use the CEs defined in Eqs. (3.2) and (3.3), i.e., Eqs. (4.94) and 
(4.95) will be assumed. Thus Eqs. (4.64) and (4.65) will also be part of the Euler a-e 
scheme. In the following, we shall describe the rest of the Eider a-e scheme. 


As a result of their definitions, evaluation of involves the inversion of the 4x4 

matrices , and , m = 1,2, 3. To simplify computation, we shall 

V J j,k \ / j,k 


assume that 


, "+ 1/2 


r j+l/3,*+l/3 


(j,Jc,n + 1/2) € ill, 


(4.96) 


mr- 

k ,+ : 

)" 

f j-2/3,k+ 1/3 ’ 

(j,k,n + 1/2) e fii, 

(4.97) 

«>r= 


)” 

J+l/3,fc — 2/3 ’ 

(j,k,n + 1/2) 6 ffy, 

(4.98) 

w- 

« 

. n+l/2 

' j-l/3,fc-l/3 ’ 

(j,k,n + 1) 6 £^ 2 ) 

(4.99) 




(4.100) 

and 

W = ’ U ’ k ' n + 1)€ ^ 

(4.101) 

By using Eqs. (4.38), (4.41), (4.44), (4.47), (4.50), (4.53), (4.74)-(4.85), and (4.96)-(4.101), 
Eqs. (4.66)-(4.71) imply that (i) 


& = « [s- (/ + f< + ) st - (i + ^ + ) < 1/w/3 . 

(4.102) 


sf ' = if' « [» + (21 -F<+) ut - (. r+F *) u-, + ]”. 2/M+1/3 . 

(4.103) 

and 

= if > « [u - (/ + F<+) ut + (2/ - Sf] J +i/3 4 _ 2/3 . 

(4.104) 
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where ( j , k,n + 1/2) G fij; and (ii) 


_ + r+ 1 / 2 

<Jy 

S' 2 <2) = jf > [a - (2/ + F< + ) tT ( + + (J - F’ + ) u+] ’ 


sT’ = »? > « [«r+ (i - f< + ) «? + (/ - f’ + ) “, + ];. 1/3 v 1/3 • 

-+] 71+1/2 

I j + 2/3,* — 1/3 


and 


5 3 (2) = 5* 3 (2) d = fu + (I - F c+ ) u + - (2/ + *”»+) u + | n+1/2 

L Jj— 1/3, 


fc+2/3 


(4.105) 

(4.106) 

(4.107) 


where (j, k, n + 1) G SI 2 . Note that sj k \ k = 1,2, and ^ = 1,2,3, respectively, are 
structually similar to s\ k \ k = 1,2, and i — 1,2,3, which were defined in Sec. 2. 

Equation (4.92) coupled with Eqs. (4.102)-(4.104) implies that 


wr-cc 

where ( j , fc,n -f 1/2) G fii, and 

(<r C* = I (**“ - 4°) • w" = | ( ? . (i) 


i+l/2 

,fe 


and (S+)jr /2 = (N° + )”t‘ /2 , (4-108) 


r(D 


) . (4.109) 


Note that: (i) As a result of Eqs. (4.102)-{4.104) and (4.109), fc (u° + )” 

can be evaluated in terms of the marching variables at the nth time level; and (ii) 
Eq. (4.109) is the Euler version of Eq. (3.20). 

Similarly, Eq. (4.93) coupled with Eqs. (4.105)-(4.107) implies that 

(“< + )ir - (^ir ’ and wc = w + )m ■ ( 4 - n °) 

where (j, k, n + 1) G a-nd 

(»r)“r=|(4 2) -4 2) ), and («r)*r = 5(4 2, -4 2) )- ( 4 - U1 ^ 

Note that: (i) As a result of Eqs. (4.105)-(4.107) and (4.111), and (u^)”^ 1 

can be evaluated in terms of the marching variables at the (n + l/2)th time level; and (ii) 
Eq. (4.111) is the Euler version of Eq. (3.36). 

Furthermore, for any (j, k, n) G ft, let (tx t )” fc denote the column matrix formed by 
m = 1,2, 3, 4. Then Eq. (4.31) coupled with Eq. (B.5) implies that 


ut = (F c+ u + + i™ + u + ) 


(4.112) 
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With the aid of Eq. (4.112), we shall carry out a development parallel to an earlier one that 
involves Eqs. (3.8)-(3.48). As will be shown, equations obtained in the new development 
are essentially the “images” of Eqs. (3.8)-(3.48) under the following mapping: 


u — ► u, u' — > u', u t — » Ut, — > F ^ + , 


/ -»/ t —t 

u c ”*■ U C’ S “*>’ 

“c -* “< + . 

< 

,,°+ _4 ,7°+ M °+ . -70+ 

U C “< ’ u v u v ’ 

“< + - J < + ' 

and 

To proceed, for any (7, Ar, n + 1/2) € Qj 

, let 



F v+ , 




u 


'+ 




(4.113) 


and 


-/n+ 1/2 def 
U i+l/3,Jfc+l/3 ~ 


— / n+l/2 def 

U j-2/3,k+l/3 “ 


-/n+1/2 def 
U /+l/3,*-2/3 “ 



n 


j+l/3,fc+l/3 


n 


/-2/3,i+l/3 


n 


j+l/3,*-2/3 


) 


7 


With the aid of Eq. (4.112), Eqs. (4.114)-(4.116) imply that 


— /n+l/2 
U j+l/3,k+l/3 


u-2 (f c+ u c + + F’ + ti+)] ” 


j+l/3,*+l/3 


(4.114) 

(4.115) 


(4.116) 


(4.117) 


and 


— / n+l/2 
U j-2/3,*+l/3 


= U- 2(F C+ U + 


+ 



n 

i-2/3,*+l/3 ’ 


(4.118) 


-/n+l/2 

U j+l/3,k-2/3 


u -2 (f c+ u + + +)] " 


j+l/3,*-2/3 


(4.119) 


For each m = 1, 2, 3,4, the earlier geometric argument involving Fig. 10(a) can be re- 
peated here with u£tj£ 2 * +1/ ,, «'"t/Y,A+i/ 3 > 311(1 vI/5-2/3’ in Fi S- 10 ( a ) bein g replaced 

by the m-th components of ^J+ 1/s , u^Jsj+i/s 311(1 “i+i/M-2/3’ respectively. This 
new argument leads to the definitions: 


wz l/2 = 


( s 


./n+l/2 

J+l/3,i+l/3 


-/n+l/2 

U j-2/3,Jk+l/3 



(4.120) 


/ \ n-4-1 /2 def 

{ Ur >)j,k 


( fi 


./n+l/2 

>+l/3,Jt+l/3 


-/n+l/2 

U i+l/3,fc-2/3 



(4.121) 
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and 


(u'+' l n+1/2 d - 

\ u i )j,k - g \ u Oi,k » 


and 


<' 7 H -' t T> + 1 / 2 def /, 7>\ n + 1 / 2 

^ Ut » t>,k g \ u v'j,k 


(4.122) 


Equations (4.120)-(4.122) are the Euler versions of Eqs. (3.16), (3.17) and (3.19), respec- 
tively. 

Next we introduce the Euler versions of Eqs. (3.21) and (3.22), i.e., 


(<*?)", t' /2 = 2 [K + )£ 1/2 - (“< + )"' 


n+1/2 

it 


, (j,Jb,n + 1/2) eft!, (4.123) 


and 


(<*OS ,/2 = 2 [« + )"t 1/2 - (u,° + )”t 1/2 ] , (;. *." + 1/2) € n,. (4.124) 


Then, with the aid of Eqs. (4.102)-(4.104), (4.109) and (4.117)-(4.122), one has 


,fc+l/3j ’ 

(4.125) 

and •• j : ,.: 


j+l/3,fc+l/3j ' 

(4.126) 

Thus both (d^^^)”^ 1/,2 and (du^)^ 1 ^ 2 are functions of the marching variables at the nth 
time level. 

Similarly, for any (j, k,n + 1) € 0,2, let 


(<«; );,r /2 = 1 (“ - 2“-< + + 4u-+) " 


>+1/3, k— 2/3 




WX ,/2 = 5 [(“+ 4 “t- 2 “t)”_ w+1/s - (“- 2S C + - 2 N + )" 


>+1/3 


■+f 

def 

+ 

At _ N 

l n+1/2 

U j-l/3,fc-l/3 


T“* y 

> — l/3,fc— 1/3 

/ n*f 1 

def 

-f 

At _ N 

n+1/2 

U j+2/3,fc-l/3 


T- y 

>+2/3, k— 1/3 

,•<( n+1 

def 

+ 

At _ ^ 

n+1/2 

U j-l/3,ik+2/3 



> — l/3,fc+2/3 


With the aid of Eq. (4.112), Eqs. (4.127)-(4.128) imply that 


*#£,»- 1/3 =[“- 2 0* + «< + + F,+ “, + )]“. 


n+1/2 


1/3, k— 1/3 


(4.127) 

(4.128) 


(4.129) 


(4.130) 
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II * H'MIBIIWIIIIIimiHllllllll 


*&}. = [ff - 2 (j=<+<r+ + J”'+<r+)] " 


n+1/2 


and 


j+2/3,*— 1/3 
n+1/2 


= [s-2(^ + \ + + ^ + “t)]”_ 1/3 , 

The Euler versions of Eqs. (3.33)-(3.35), (3.37), and (3.38) are 

/'i 7 , ') n + 1 ^ /, 7 'n+l 7 / n+l \ /a/- 

( U C'j,* — (^ U i+ 2/3, *-1/3 U i — 1/3,* — 1/3 J / A C’ 


*+2/3 


<'iT'l n+1 — f (u ,n+1 -u' n + 1 ^ / A n 

~ — 1/3,*+2/3 U j-l/3,k-l/3j / Ar ?’ 


(4.131) 

(4.132) 


(4.133) 

(4.134) 


(ffDiJ 1 = x( s <)m’’ “ d W “ xW** 1 ’ (4135) 


and 


(<*?)£' H ? 2 [(«'+)”+■ - . 0\M + 1 ) € n 2 , 

(<*<)?,:* = 2 [(ff; + )”,r - (ff; + )”t‘] , 0, t,«+i)e « 2 . 


Combining Eqs. (4.105)-(4.107), and (4.130)-(4.137), one has 

GW = 5 [(“ + 2 “t + 2 <)"-! / M- 1/ 3 - + 2ff » + )," 

and 

^ 1 T/ ,\n+l/2 / . ,\ n+1/2 

W)?* =0 ( « + 2u^ + + 2u + ) - ( « + 2u + - 4u + ) 

V V ' J ’ k 3 V < * Jj-l/3,k-l/3 \ c 1 / 3 , 


n+1/2 

j+2/3,*— l/3j 


*+2/3 J 


(4.136) 

(4.137) 


, (4.138) 


. (4.139) 


Thus both (du£ )”|‘ 1 and (du+ )"£ 1 sure functions of the marching variables at the (n+l/2)th 
time level. 

The Euler a-e scheme can now be stated using the above definitions. It consists of 
two marching steps. The first is formed by Eq. (4.64), 


and 



n+1/2 

n+1/2 

j.* 



n+1/2 

j,k 


= (“7 + );,r /s 



n+1/2 
j,k » 


+ e(du$) 


n+1/2 

J,* 


1 


(4.140) 

(4.141) 
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where (j,k,n + 1/2) 6 fti, and e is an adjustable parameter. It was explained earlier 
that the expressions on the right sides of Eqs. (4.140) and (4.141) are functions of the 
marching variables at the nth time level. Moreover, according to Eqs. (4.123) and (4.124), 
Eqs. (4.140) and (4.141) can also be expressed as 



« i/2 = ( s < + >s i/2 + ( £ - 1/2 

(4.142) 

and 

w )"?' 1 = wit 1 '’ + (* - i/ 2 )(<j<)“,: ,/2 , 

(4.143) 

respectively. 

The second marching step is formed by Eq. (4.65), 




(4.144) 

and 


(4.145) 


where (j, k, n + 1) £ f& 2 - It was explained earlier that the expressions on the right sides of 
Eqs. (4.144) and (4.145) are functions of the marching variables at the (n + l/2)th time 
level. Furthermore, according to Eqs. (4.136) and (4.137), Eqs. (4.144) and (4.145) can 
also be expressed as 


and 




\n+l 

)j,i 


n+l 


+ (e-l/2 )(*?)£ 


(4.146) 



n+l 

i.* 




(4.147) 


respectively. 

Note that, because of the assumptions made in Eqs. ( 4.96)-(4.101 ), the Euler a scheme 
is not the special case of the Euler a-e scheme with e = 0. 

Finally we shall construct the weighted-average Euler a-e scheme, i.e., the Euler version 
of the weighted-average a-e scheme. The development follows a line of argument similar 
to that used in the construction of the weighted- average a-e scheme. Thus only the key 
definitions will be given. 

For any ( j,k,n + 1/2) € fii, the Euler versions of Eqs. (3.56)-(3.66), (3.69)-(3.74), 
and (3.80) are 


del jln+1/2 .-rn+l/2 

Xl ~ U j+l/3,k+l/3 U j,k > 


(4.148) 


def 

x 2 — Uj_2/3,Jk+l/3 


, V n + 1 /2 

J.* ’ 


U 


(4.149) 
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_* dcf -*/ n-f 1/2 -*n+l/2 

X 3 ~ U i+l/3,fc-2/3 — U i,ifc ’ 

(4.150) 


(ff< (,, )S ,/, = -(2*« + *»)/*<. 

(4.151) 


( ff » 1) )J» 1/s =* _ (ft + 2ft)/ai, 

(4.152) 


(“c 2) )7,r /2 = (2^ + 

(4.153) 


(sf);,t 1/2 = (ft-ft)/iv, 

(4.154) 


(“c 3> )”,i 1/2 = (ft - ft)/*c> 

(4.155) 


(“1 3) )"t 1/2 = (2ft + ft)/A^. 

(4.156) 


(“1 <> )m ,/2 I + (“T)m I/2 * = 1.2,3, 

(4.157) 


(“f ) );,r /2 |+(“T);,r /2 |. < = 1 . 2 , 3 , 

(4.158) 


(ft <1> ) J T /2 = -£;(ft + ft). 

(4.159) 


/Wi>\r+i/2 (36 + u>)x 2 + (36 - u>)x 3 

l “» 'l'* “ 2t»ft 

(4.160) 


w /2 = §, 

(4.161) 


/^(2)x«+i/2 (36 + ty)£i + 2 iux 3 

(4.162) 


^ h* 2w ’ 

(4.163) 


t -(3^n+l /2 (™ ~ 36)xi + 2WX 2 

{U » '** ~ 2uifc 

(4.164) 
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and 


(< 0+ )m 1/2 = and (ff<« + );,t I/2 = )"t 1/2 , (4.165) 

respectively. 

Let the mth components of the 4x1 column matrices (u^ + )""£ 1 ^ 2 , (u I ^ + )J fc 1 ^ 2 , 
(-W ) n+1 /2 J md („-<<>)»+■/*, be denoted by (u£ + )£ 1/2 , (ni£, + )”t’ /2 . and 

(umj);^ 2 , respectively. Then, for m = 1,2,3, 4, and i = 1,2,3, the Euler versions of 
Eqs. (3.79), (3.85), and (3.86) are 


K, = v/(»&) 1 + («!?,)», 


(4.166) 


0, 

U ™< = ) (gmi*g3 )° «S + +(«»»«»! )“ ttj% + + (» ml »«,)- »g + 


(^ml^m2)“ + ($m2$m3 )° + (^m3^ml) a 


0, 


if $ml — ^m2 — $m3 — 0> 

, otherwise, 

(4.167) 

if #ml = ^m2 = ^m3 = 0; 

. . . v . . . « . . _ , , otherwise, 

(0ml0m2) a + (^m2^m3)“ + (^m3^ml)“ 

(4.168) 

respectively. Here (i) a > 0, and (ii) each symbol in the last two equations is associated 
with the mesh point (j,k,n + 1/2). 

Next we consider the case with (j, k, n + 1) € Q 2 - The Euler versions of Eqs. (3.90)- 
(3.106) sire 

r/3-^r, (4.169) 


and 


tc-J- def 

u = < 

TOf/ ' 


(^m2^mj)° Um^ + +(^m3^ml) a +(^wl^m2)° 


,7_ ^ 

- U >+2/3,fc-l/3 U J,fc ’ 


(4.170) 


,7, ^ it' n + 1 _ T 7 n + 1 

— U ;-l/3,fc+2/3 U j,k > 


(ffW^+l HF (2yj + jr,)/iC. 


(4.171) 

(4.172) 
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and 


(N (1) ®‘ d = (?2 +2J 3 )/Ar(, 


(\ (2) );,t ld = -(2w+fi)/A(, 


(“J 2 ) )",t‘ =' (ft - 

(“c 3) )?.t‘ “(ft-ft)Mc, 

')#’ = -(2ft + ft)/**. 


(„«)"+■ 4* («r‘0)j+* £ + £, 


a< 


^7 




9x 


t = 1,2,3, 


(W*(«< M «'g+(W& 


ay’ 


^ = 1,2,3, 




(4.173) 

(4.174) 

(4.175) 

(4.176) 

(4.177) 

(4.178) 

(4.179) 

(4.180) 


,-riK«+i ( 36 + W )V 2 + (36 - *»)& 
>i* ~ 2 wh 


(4.181) 


(*i 2) )Jf = 


3yi 
2w ’ 


(4.182) 


,.-W2Kn+l ( 36 + w )Vl + 2w & 

{U * >i’ k ~ 2tvh 


(4.183) 


«w - -i- 


and 


(w - 3fe)yi + 2wy 2 


= 2w/r 


(4.184) 


(4.185) 


respectively. Moreover, with the understanding that all symbols are associated with the 
mesh point (j, k,n + 1) € O 2 , Eqs. (4.165)-(4.168) remain valid for the current case. 
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The weighted-average Euler a-e scheme can now be stated. It consists of two marching 
steps. The first is formed by Eq. (4.64), 



( s < + )m’ / 2 = (“Di: i/2 + (« - i/2x«w?)”t ,/2 . 

(4.186) 

1 

and 

W)”t' /2 = («“ + )”+ 1/2 + (e - l/2)(iu+)”t ,/2 , 

(4.187) 


where (j, k,n +• 1/2) € ftj. The second is formed by Eq. (4.65), 



(u+)”^ = (S { "+)J+‘ + (« - 1/2)(<H+)S , ) 

(4.188) 

: 

and 

K + )"t‘ = + (« - 

(4.189) 

- 

where (j, fc, n + 1) € 

Because (i) a fractional power is costly to evaluate, and (ii) evaluation of (6 m e) a does 
not involve a fractional power if a is an even integer, the weighted-average Euler a-e scheme 
is more computationally efficient if a is an even integer. 
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5. Stability Analysis 


The stability of the a, the a-/x, and the a-e schemes will be studied using the von 
Neumann analysis. For all (j, k, n ) € Cl, let 

q(j,k,n) = q*(n,9 <: ,6 v )e ,b0 < +k *'\ (t d =v^l, -ir < 6{,d n < ir). (5.1) 

where q*(n,9^,9 v ) is a 3 x 1 column matrix. Substituting Eq. (5.1) into Eqs. (2.84) and 
(2.85), one obtains 

f(n + m — 1/2, 6 ( , «,) = 0,)A/< 2 >(0 c , «,)] "V(n - 1/2, 0 C , *,), (5.2) 

and 

r(» + m,« c ,«,)= «,)«<«(«(, «,)] m f (n, «<,«,), (5.3) 

respectively. Here (i) n = 0, ±1, ±2, . . .; (ii) m = 0, 1,2,.. .; and (iii) 

M (1) (0 f ,0„) d =l f + ga) c («73)(-2« <+ ^) + Q(i) e «/3)(*<-2*,), (5.4) 

and 

M (2) (0 C ,0,,) d = g5 2) e _(,/3)( ® <+ *’’ ) + Q 5 2) C (,/3)( “ 2tf<+d,) + g( 2) e -(*'/ 3) (^- 2 ^). (5.5) 

Note that Eqs. (2.84) and (2.85) are valid for the above three schemes if Q^\ k = 1,2, 
and i — 1,2,3, are defined using (i) Eq. (2.80) for the a scheme, (ii) Eq. (2.107) for the 
a-jj. scheme, and (iii) Eqs. (3.50)-(3.55) for the a-e scheme. Equation (5.2) implies that 
the amplification matrix among the half-integer time levels is 9 r) )M^ 1 \9^, #,); 

while Eq. (5.3) implies that the amplification matrix among the whole-integer time levels 
is 

According to a theorem given in Appendix C, the above two amplification matrices 
have the same eigenvalues. These eigenvalues may be referred to as the amplification 
factors. The amplification factors axe functions of phase angles 9^ and 9 V . In addition, 
they are functions of a set of coefficients which axe dependent on the physical properties 
and the mesh parameters. These coefficients axe (i) and u v for the a scheme; (ii) v^, 
z/,, and e for the a-e scheme; and (iii) z^, u v , £ v , and £ r for the a-fi scheme. Let A 1? 
A 2 , and A 3 denote the amplification factors. In the present paper, a scheme is said to be 
stable in a domain of the above coefficients if, for all coefficients belonging to this domain, 
and all 6$ and 9 V with —ir < 9(,9 V < ir, 

|Ai|<l, |A 2 | < 1, and |A 3 | < 1. (5.6) 

Consider the a scheme. By using its two-way marching nature (see Fig. 12), it is 
shown in Appendix C that, for any given z^, v v , 9^, and 9 V , 

|A 1 A 2 A 3 | = 1. (5.7) 
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It follows from Eqs. (5.6) and (5.7) that the a scheme must be neutrally stable, i.e., 

I Ai| = JA 2 1 = | A 3 [ = 1, ^ ^ —ft < 0c 0 n <tt, (5.8) 

if it is stable. In other words, the a scheme is free of numerical diffusion [5, p.18] if it is 
stable. Moreover, a systematic numerical evaluation of Ai, A 2 , and A 3 , for different values 
of V{, v v , 8(, and 0 V , has confirmed that the a scheme is indeed neutrally stable in the 
stability domain defined by Eq. (2.110). In the following, we shall discuss the meaning of 
this stability domain. 

Let ( j , k,n + 1/2) 6 Cl\. According to Eq. (2.84), the marching variables at (j, k,n + 
1/2) are completely determined by those of seven mesh points at the (n — l/2)th time 
level. According to Fig. 13(a), one of them, i.e., the mesh point (j, k, n — 1/2), is located 
directly “below” the mesh point (j,k,n + 1 /2). The other six are the vertices of a hexagon. 
As a result, in this paper, the interior and boundary of the hexagon shall be considered as 
the numerical domain of dependence of ( j , fc, n + 1/2) at the ( n — l/2)th time level. Note 
that the dashed lines depicted in Figs. 13(a) and 13(b) are the spatial projections of the 
boundaries of CEs. 

The a scheme is designed to solve Eq. (3.1). For Eq. (3.1), the value of u is a constant 
along a characteristic line. The characteristic line passing through the mesh point (j, k,n+ 
1/2) will intersect a point on the plane t = t n-1 / 2 . The latter, referred to as the backward 
characteristic projection of the former at the (n — 1/2) th time level, is the “domain” of 
dependence of the former at the (n — l/2)th time level. It is shown in Appendix C that the 
backward characteristic projection is in the interior of the numerical domain of dependence 
if and only if Eq. (2.110) is satisfied. 

Let (j, k,n + 1) € and consider Fig. 13(b). Using a line of argument similar to 
that presented above, it can be shown that the backward characteristic projection of the 
mesh point (j, k,n + 1) at the nth time level is in the interior of the numerical domain of 
dependence at the nth time level if and only if Eq. (2.110) is satisfied. 

At this juncture, note that the concept of characteristics was never used in the design 
of the a scheme. Nevertheless, its stability condition is completely consistent with the 
general requirement that an explicit scheme for Eq. (3.1) is stable when the domain of 
dependence of Eq. (3.1) is a subset of the numerical domain of dependence. 

Next we consider the stability of the a-e scheme. Recall that the 1-D a-e scheme [5] is 
not stable for any Courant number v if e <0, or e > 1. Similarly, the results of numerical 
experiments indicate that the current a-e scheme, except for some possible isolated points, 
is not stable in any domain on the v^-v^ plane if e < 0 or e > 1. For any e with 0 < e < 1, 
the a-e scheme has a stability domain on the plane. The stabilibity domains for 

several values of e were obtained numerically. As shown in Figs. 14(a)-(c), these domains 
(shaded areas) vary only slightly in shape and size from that depicted in Fig. 9. They 
become smaller in size as e increases. ~ 

Let Ai, A 2 , and A 3 be defined such that 

|A 3 | < |A 2 | < |A a |. (5.9) 
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Then Ai can be referred to as the principal amplification factor; while A2 and A3 referred 
to as the spurious amplification factors [1]. In general, the principal amplification factor is 
the deciding factor in determining the accuracy of computation [1]. Particularly, numerical 
solutions may suffer annihilations of sharply different degrees at different locations and 
different frequencies if numerical diffusion associated with Aj varies greatly with respect 
to 9(, 6 V , and u n [5, p.20]. Assuming Eq. (5.6), then (1 — |A*|) is a measure of the 
numerical diffusion associated with A*, £ = 1,2,3 [5, p.18]. For a given e, let D(e) denote 
the stability domain of the a-e scheme on the u^-i plane. Let 


X/(e) = f max (1 - |A*|), 

("o* ',)eD(0 


£ = 1,2,3; 0 < e < 1. 


(5.10) 


Then, for a given e and each £, (1 — |A/|) is bounded uniformly from above by Xt( 6 )- The 
numerically estimated values of x*( e ) 3X6 plotted in Fig. 15. From this figure, one con- 
cludes that the numerical diffusion, particularly that associated with Ai , can be bounded 
uniformly from above by an arbitrary small number by choosing an c small enough. Note 
that this property is also shared by the 1-D a-e scheme (see Eq. (3.19) in [5]). Moreover, 
the results shown in Fig. 15 indicate that X 2 ( e ) and X 3 ( 6 ) 3X6 much larger than Xi( e ) m 
the range of 0 < e < 0.5. Thus, in this range, the spurious part of a numerical solution is 
annihilated much faster than the principal part. Also it is seen that the numerical diffusion 
associated with the principal solution, measured by Xi( e )s increases with e in the range of 
0 < e < 0.7. 

Finally we discuss the stability of the a-fi scheme. Note that this scheme is not defined 
if A^ = 0 or = 0. One form of A^ 1 ^ and A^ is given in Eqs. (2.86) and (2.87). Another 
form is given in Eqs. (A. 15) and (A. 16). Let fi — 0. Then A^ = 0 or a^ = 0 occurs 
only on the six straight lines on the 1/^-1/, plane which are depicted in Fig. 16. The shaded 
area depicted in the same figure is the region that satisfies Eq. (2.88). It was explained in 
Sec. 2 that the curves of singularity on which A^ = 0 or A^ = 0 cannot enter the shaded 
region if n > 0. In general, for a given set of £< > 0, £, > 0 and £ r > 0, the v^-v v plane 
can be divided into two regions by the curves of singularity. The “inside” region R\ is the 
maximal connected open set on the plane that contains (i) the shaded area depicted 
in Fig. 16, and (ii) no point at which A^ 1 ' = 0, or A^ = 0. The “outside” region R 2 is the 
rest of the v^-u n plane. 

For the a-p scheme, A/, £ = 1,2,3, are functions of the phase angles 0^ and 9 V , and 
the coefficients £,, and £ r . For a given set of and the stability domain 

on the plane generally covers part of R\ and part of i?2- In the following discussion, 
only the stability domains in i?i are considered. 

For the special case with a£ = A 7 ? = at, (<;, £ 7 , and £ r share a common value, say £. 
The stability domains (shaded areas) for £ = 10“ 5 , £ = 10~ 3 , and £ — 0.1 are plotted in 
Figs. 17(a)-(c), respectively. 

Next we assume that a( = A 77. Then Let a be the angle formed by the sides 

DB and DF which are depicted in Figs. 7(a)-(c). Then Eq. (2.32) implies that 


fr = 2(1 — cosa)£, 


(5.11) 
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where £ is the common value of ^ and The stability domains for two pairs of (£,<*) 
are depicted in Figs. 18(a) and (b). • • ■ . .. _ 

From the results shown in Figs. 17 and 18, and the results of other numerical exper- 
iments, it appears that the a-/x scheme is unconditionally stable when = u v — 0. The 
last condition is equivalent to = a n = 0 or a x = a y = 0. 
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6. Numerical Results 


In [8], several numerical solutions of Eqs. (2.1) and (3.1) generated using the a-p and 
the a-e schemes are compared with the exact solutions or the numerical solutions generated 
using traditional methods. These comparisons show that the a-e scheme, which includes 
the a scheme as a special case with e = 0, is an accurate solver for Eq. (3.1). They also 
show that the a-fi scheme can obtain highly accurate solutions of Eq. (2.1) as long as 
the viscosity coefficient p is not too large. Note that a convection-diffusion probelm is 
fundamentally an initial- value/boimdary- value problem. The current explicit a-fi scheme 
obviously cannot model such a problem unless the contribution of the diffusion term is 
small compared to that of the convection term. 

The a-e scheme was also generalized in [8] to solve the 2-D inviscid Burgers’ equation. 
In spite of its simplicity, particularly the fact that it does not use (i) any mesh refinement 
technique, and (ii) any moving meshes, this new solver is capable of generating highly 
accurate shock solutions. The shock discontinuities are alm ost, resolved within one mesh 
intervals. 

In this section, accuracy of the weighted-average Euler a-e scheme defined in Sec. 4 will 
be evaluated using a steady-state shock reflection problem [24]. The computation domain 
and the shock locations ( AE and EF) are depicted in Fig. 19. The lower boundary is a 
solid wall. Assuming 7 = 1.4, the exact Euler solution to this problem is: 

(a) In the region ABE , 


u = 2.9, 

v = 0., p — 1.0, p = 1. 0/1.4. 

(6.1) 

(b) In the region AEFD, 



u = 2.6193, v 

= -0.50632, p = 1.7000, p = 1.5282. 

(6.2) 

(c) In the region ECF , 



u = 2.4015, 

v = 0., p = 2.6872, p = 2.9340. 

(6.3) 


Note that the Mach number is equal to (i) 2.9 in the region ABE ; (ii) 2.3781 in the region 
AEFD ; and (iii) 1.9424 in the region ECF. 

The mesh used in the current numerical calculations is depicted in Fig. 20. Again 
a mesh point £ fii is marked by a solid circle; while a mesh point £ f) 2 is marked by 
an open circle. The mesh is a special case of that depicted in Figs. 1-4 with 6 = 0. 
Note that (i) only the mesh points £ are present at the inflow boundary, and (ii) the 
mesh parameter w is so chosen that only the mesh points € ft 2 are present at the outflow 
boundary. Moreover, for simplicity, a mesh point and the corresponding marching variable 
will be identified by the time-level number n, and two new mesh indices r and s which are 
given in Fig. 20 as a pair of integers enclosed in a parenthesis. Note that, for the mesh 
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points E fii , r = 1, 2, 3, . . . , R, R + 1, and s = 1, 2, 3, . . . , S. On the other hand, for the 
mesh points € CI 2 , r = 1, 2, 3, . . . , R, R + 1, and s = 1, 2, 3, . . . , S, S + 1. Obviously two 
different mesh points at the same time level always have different pairs of r and s. 

In the current numerical calculation, at the time level n = 0, u m , m = 1, 2, 3,4, at all 
mesh points are calculated using Eq. (6.1). Also we assume that 


u rn< = “«.,= °» m = 1, 2, 3, 4. (6.4) 

The above initial conditions are also assumed at the inflow boundary for all n = 1, 2, 3, — 

At the upper boundary, for all n = 1/2, 1, 3/2,2, . . ., Eq. (6.4) is also assumed. Moreover, 

u m , m = 1,2, 3, 4, are calculated using Eq. (6.2). t “ ' j 

To impose the proper b ound ary conditions at the lower boundary, note that the solid j 

wall boundary conditions at BC (see Fig. 19) are equivalent to the condition that the flow 
field below BC is the mirror image of that above. By using Eq. (4.1) and the fact that 
y = 0 at any point on BC, it can be shown that the last condition implies that 

«m(*, ~y) = u m (x,y), m = 1,2,4, and u 3 (x, -y) = -u 3 (x,y). (6.5) 


and 


du m (x,-y ) du m (x, y) ^ du m (x,-y ) du m (x,y ) __ , , 

= , and — = — , m = 1,2,4, 

Ox ox oy oy 

du 3 (x,-y) du 3 (x,y ) du 3 (x,~y ) _ du 3 (x,y) 

dx dx 5 311 dy dy 


( 6 . 6 ) 

(6.7) 


Consider the mesh depicted in Fig. 20. Then it becomes clear that the numerical analogues 
of Eqs. (6.5)-(6.7) are 


( U ")h+M = Mr, s' m=1 ’ 2 ’ 4 ’ and (“3 )r +M = -(«3)r )S , 


( 6 . 8 ) 


( U "»*)r+ 1,3 “ ( U mx) ( u my) R+1( , — ( u my) Ra , m— 1,2,4, (6.9) 

and 

( U 3*)R+M=-( W 3*)H, a > and ( U 3 »)r + 1) , = («3,) R)i , (6.10) 


respectively. According to Fig. 20, the range of s in Eqs. (6.8)-(6.10) is dependent on 
the time level n. Let (i) S + S + 1, and S~ *== S if S is even; and (ii) S + d = S, 

and S~ d = S — 1 if S is odd. Then (i) s = 2, 4, 6, . . . , S~ if n = 1/2, 3/2, ..., and (ii) 
s = 1, 3, 5, . . . , 5 + if n = 1,2, — Furthermore, by using Eq. (B.4) with b = 0, it can be 
shown that Eqs. (6.9) and (6.10) are eqivalent to 



n 

R+h* 


= to.,)*., - 


and 





m = 1,2,4, 


( 6 . 11 ) 
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ii min uni 1 1 


and 


( 6 . 12 ) 



n 

J£-f l,a 




“ d K)«+l,.= 



5 


respectively. Equations (6.8), (6.11) and (6.12) axe the boundary conditions at the lower 
wall (a solid wall) in the current numerical calculations. In other words, the marching 
variables associated with the mesh points below the solid wall will be determined using 
Eqs. (6.8), (6.11), and (6.12). 

Next we discuss the outflow boundary conditions. For any n = 1,2,3,..., and r = 
1 , 2, 3, . . . , R, we assume that 


(«»)" s+ , = («m)" s I/2 , m = 1,2, 3, 4, (6.13) 

(“">*)"s+l = 0. m ~ 1,2, 3, 4, (6.14) 

and 

(“»»)".*+ 1 = (“-»)"s 1/2 . m = 1,2, 3, 4. (6.15) 

When the time-marching solution reaches its steady-state limit, the above conditions can 
be considered as a result of the requirement that the partial derivatives of the flow variables 
with respect to x are zero at the outflow boundary. By using Eq. (B.4) with b = 0, it can 
be shown that Eqs. (6.14) and (6.15) are equivalent to 

= m = 1 ’ 2 ’ 3 ’ 4 ’ < 616 ) 

and 

(«m„)r,S+l = \ (“m, “ ’ TO = 1 ’ 2 > 3 ’ 4 ’ ( 6 - 17 ) 

where n = 1,2,3,..., and r = 1,2,3,...,#. Equations (6.13), (6.16), and (6.17) axe 
the outflow boimdaxy conditions in the current numerical calculations. As a result, the 
marching variables at the outflow boundary will be determined using these equations. 

With the add of the above initial and boundary conditions, the marching variables at 
all time levels can be determined using the weighted-averaged Euler a-e scheme. As an 
example, at any n = 1/2, 3/2, . . ., the marching variables associated with the mesh point 
(2, 1) (marked by a solid circle in Fig. 20) can be determined in terms of those associated 
with the mesh points (1, 1), (2, 1), and (2, 2) at the (n — l/2)th time level (maxked by open 
circles). As another example, at any n = 1,2,3,..., the marching variables associated 
with the mesh point (1,3) (marked by an open circle) can be determined in terms of 
those associated with the mesh points (1, 2), (2, 3), and (1,3) at the (n — l/2)th time level 
(marked by solid circles). 

According to Fig. 19, the distance between the inflow and the outflow boundaries is 
4., while the distance between the upper and the lower boundaries is 1.. On the other 
hand, according to Fig. 20, the above two distances are w ■ S and 2 h ■ R, respectively. Thus 

w = and h = 2 ( 6 - 18 ) 
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Because 6 = 0, the geometric parameters w, h, and 6 are determined if R and 5 are given. 
In addition to the initial conditions, the boundary conditions, and the integers R and 
S, the other input parameters for the current numerical calculations are e, a, At, and a 
positive integer n t . Here we assume that the time marching ends at the n t th time level, 

i.e., at t = T A = n t - At. 

It is shown in Appendix D that, for any Euler solver constructed in Sec. 4, a local 
CFL number v t associated with any mesh point (j, &, n) € ft can be defined in terms of 
u, v, c, w, h and At. Here u, v, and c are the x-velocity, the y-velocity, and the sonic 
speed at the mesh point, respectively. Two global CFL numbers are considered in the 
current calculations. The first, denoted by v ema , is the maximum of v t with respect to the 
steady-state solution given in Eqs. (6.1)-(6.3). The second, denoted by u em , is the largest 
value of v e ever reached at any mesh point (j, k, n) € ft, where n = 0 , 1 / 2 , 2 , 3/2, 2 , . . . ,n*. 
Excluding the initial and the boundary conditions, v em $ i s dependent on R, S, and At 
only. On the other hand, u tm is a function of R, 5, At, n t , e, and a. According to a series 
of numerical experiments, the value of a plays only a minor role on the stability of the 
weighted-average Euler a-e scheme. Generally the scheme is stable if 

0 < e < 1 and v tm < 1 . (6.19) 

To measure the convergence of a time-marching solution to the corresponding steady- 
state solution (note: this steady-state solution generally differs from the exact solution 
given in Eqs. (6.1)-(6.3)), for any n = 1 , 2, 3, . . . , n t , and m = 1,2, 3, 4, let 



f 1 

’s+i 

r(i)+R-l 

1 

log 10 < 

1 1 
I RSCm 

£ 

E IK.)?,. 1 



1 Nl 

*=2 

r=r(a) 

4 


( 6 . 20 ) 


Here, for any m = 1,2, 3,4, c m is the maximal value of |u m | within the exact steady-state 
solution defined by Eqs. (6.1)-(6.3). It can be shown that c\ = 2.6872, c 2 = 6.4534, 
C 3 = 0.86073, and C 4 = 15.084. Moreover, 



if s is odd; 
otherwise. 


( 6 . 21 ) 


According to Fig. 20, the summation which takes place in Eq. (6.20) involves all the mesh 
points at the nth time level excluding those located (i) at the inflow boundary, (ii) at 
the upper boundary, and (iii) below the lower boundary. Because n = 1,2,3, ... ,n t , the 
mesh points involved in summation are all marked by open circles in Fig. 20. The values 
of u m at the inflow and the upper boundaries do not change with time, while those at 
the mesh points in (iii) are dependent on the values of u m at other interior mesh points. 
Note that the values of u m at the outflow boundary change with time and are dependent 
on those at a lower time level. Because the summation involves a total of R x S' mesh 
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points, the result of this summation divided by R x S is the average value of the change 
of u m at the same spatial mesh point (measured by the absolute value of this change) 
from the (n — l)th time level to the nth time level per mesh point. This average value is 
further normalized using the constant c m . If we further assume that the time-marching 
solution converges to a steady-state solution which is similar to the exact steady-state 
solution (such that the normalization by c m makes sense), then E m (n) more or less can 
be interpreted as the average number of correct significant figures in u m at the nth time 
level as compared with the converged value of u m (which, of course, is not identical to that 
given in Eqs. (6.1)-(6.3)). 

Because the time marching solution can not reach a steady-state solution before the 
boundary conditions are fully felt at all interior points, rapid convergence generally can not 
occur before the time has elapsed which allows a fluid particle to travel the full length of 
the computation domain. It can be shown that, for the solution given in Eqs. (6.1)-(6.3), 
the average value of u over the computational domain is 2.6261. Thus an average fluid 
particle requires 4.0/2.6261 = 1.5232 time units to travel from the inflow boundary to the 
outflow boundary. The number of time steps corresponding to the above number of time 
uints is 


n c 


def 1.5232 

= At ' 


( 6 . 22 ) 


i.e., rapid convergence can not occur before n > n c . 

With the above preliminaries, the numerical results generated using the weighted- 
average Euler a-e scheme can now be presented. Six test problems, with different combi- 
nations of e, a, R, S, At, and n t , are defined in Table 1. For each problem, the values 
of T, i/ em „ v cm , and n c are also given in the same table. In Figs. 22-27, the numerical 
results (triangular symbols) of the pressure coefficient c p at n = n t for Problems #l-#6 
are compared with the exact solution (solid lines). Here 


djsf 2 f p \ 

Cp l M lo \P°o ) ’ 


(6.23) 


with Mqo = 2-9 and poo = 1.0/1. 4 being the inflow Mach number and pressure, respectively. 
Note that: (i) at the mid-section of the computation domain (y = 0.5 in Fig. 19), two 
neighboring mesh points at the same time level are separated by a distance = 2 w, and (ii) 
the mesh points at the n«th time level are marked by open circles in Fig. 20 because n t 
is a whole number. In Figs. 22-27, the values of E m (n), m = 1,2, 3, 4, are also plotted 
against n for all six test problems. In Fig. 28, twenty-six pressure contour levels between 
the values of 0.6 and 3.1 with uniform increment 0.1 were used for the contour plots of 
Problem #3. Finally, for Problem #3, a 3-D pressure-distribution plot is shown in Fig. 29. 

The significance of the results shown in Figs. 22-29 is discussed in the following 
remarks: 

(a) From Table 1 and the results shown in Figs. 22-24, it appears that the convergence to 
steady-state is much faster with a smaller value of i/ em (or v ems ). As a matter of fact, 
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convergence to steady-state can reach a plateau representing some number of correct 
significant figures if v em is too close to 1. Prom Table 1 and a comparison among 
Figs. 22, 25, and 26, one also concludes that slower convergence generally occurs with 
a value of e much smaller than 0.5. A comparison between Figs. 22 and 27 reveals 
that a change of the value of a from 2 to 1 also cause a slight decrease in convergence 
rate. Because numerical diffusion generally increases with (i) a smaller value of v tm , 
(ii) a larger value of e, and (iii) a larger value of a, one may conclude that faster 
convergence generally occurs with larger numerical diffusion. This trend is consistent 
with the fact that shocks cannot be formed without physical or numerical diffusion. 

(b) The effectiveness of the weighted- averaging as a tool to surpress numerical oscilla- 
tions near discontinuities is clearly demonstrated by the results shown in Figs. 22-29; 
Moreover, the present weighted-averaging does not cause the smearing of shock dis- 
continuities and has no discernible effect on the smooth part of the solution. From 
table 1 and a comparison between Fig. 22 and 27, one also concludes that the increase 
of the value of a from 1 to 2 has a marginal impact on the numerical result?. 

(c) Comparing the numerical results shown in Figs. 22-27 with the exact solution, one 
concludes that the weighted-average Euler a-e scheme is capable of generating highly 
accurate solutions for the steady-state shock reflection problem under consideration. 
Also a comparison of the results shown in Figs. 22, 23, and 25-27 reveals that accuracy 
of the numerical results generally is not sensitive to the change of the values of i/ em , 
£, and e. An exception is that numerical results may become more diffusive and thus 
shock resolution becomes less sharp if the value of e is too large, e.g., e = 0.8 in 
Problem #5. Finally, a comparison of the results shown in Fig. 24 (Problem #3) with 
the results of other test problems reveals that accuracy increases sharply with the 
decrease of the mesh size. It is seen that both the primary and the reflected shocks 
are resolved by a single data point in Fig. 24. 
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7. Conclusions and Discussions 


A new numerical method is being developed for solving one- dimensional and multidi- 
mensional flow problems. This new method represents a clear break from the traditional 
methods in the basic concept of numerical discretization. It emphasizes simplicity, gener- 
ality, and accuracy. The history of this new method and the considerations that motivates 
its development are clearly described in Sec. 1. 

In this report, we explain how the same set of design principles which was used to 
construct several solvers for 1-D time-marching problems [5] can be used to construct 
their 2-D counterparts. Because of the similarity in their design, each of the present 2- 
D solvers virtually shares with its 1-D counterpart the same fundamental characteristics. 
Furthermore, it has been shown that the 2-D solvers, as in the case of the 1-D solvers, 
generally are more accurate than the traditional solvers in spite of the advantage the 
present solvers have over the latter in simplicity and generality. Accuracy of the present 
2-D Euler solver is most vividly demonstrated by the pressure- contour plot (Fig. 28) and 
the 3-D pressure-distribution plot (Fig. 29) it generates for a famous shock reflection 
problem [24]. Both the primary and the reflected shocks are resolved by a single data 
point without the presence of numerical oscillations near the discontinuity. 

The construction of the 1-D solvers referred to above is simplified by the use of a 
mesh which is staggered in time [5]. Its use results in the simplest stencil possible, i.e., a 
triangle in a 2-D space-time with a vertex at the upper time level and other two at the 
lower time level. Similarly, the construction of the present 2-D solvers is simplified by the 
use of a nontraditional space-time mesh which is also staggered in time (Figs. 1-4). Its 
use results in the simplest stencil possible, i.e., a tetrahedron (Fig. 8) in a 3-D space-time 
with a vertex at the upper time level and other three at the lower time levels. 

The meshes used by the 1-D and the 2-D solvers consist of whole-integer and half- 
integer time levels with a half-integer time level being sandwiched between two whole- 
integer time levels and vice versa. The spatial positions of the mesh points at a whole- 
integer (half-integer) time level coincide with those at another whole-integer (half-integer) 
time level. However, the spatial positions of the mesh points at a whole-integer time level 
shift from those at a half-integer time level. For the mesh used by the 1-D solvers, the 
spatial projection of a mesh point at a whole-integer time level is right at the center of 
those of two neighboring mesh points at a half-integer time level and vice versa [5]. It 
follows that the stencil of the 1-D solvers is always an isosceles triangle, i.e., one cannot 
distinguish a stencil with its upper vertex at a whole-integer time level from another with 
its upper vertex at a half-integer time level. As a result, each of the 1-D solvers constructed 
in [5] is formed by two identical marching steps. Contrarily, for the present 2-D solvers, a 
stencil (a tetrahedron) with it vertex at a whole- integer time level is different from another 
with its vertex at a half-integer time level (Fig. 8). Thus each of the present 2-D solvers is 
formed by two distinctly different marching steps. In spite of their structural differences, 
the last two marching steps compensate each other and its combination results in severed 
important symmetric properties which were discussed in Sec. 5 and Appendix C. 
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The Euler a scheme constructed in Sec. 4 is free from numerical diffusion when it 
is stable. This scheme is a limiting case of a Navier-Stokes solver currently under devel- 
opment, i.e., the former is a special case of the latter when the viscosity vanishes. As 
a result, the new Navier-Stokes solver will have a special property that a classical solver 
lacks, i.e., as the physical diffusion (viscosity) approaches zero, so does the numerical dif- 
fusion. Without this property, n umerical dissipation may overwhelm physical dissipation 
and cause a complete distortion of solutions for problems with small viscosity. Because 
a Navier-Stokes problem fundamentally is an initial- value /boundary- value problem, i.e., 
information from any spatial point can be felt instantly by other spatial points, the new 
Navier-Stokes solver is implicit when the viscosity is present. However, it is reduced to the 
Euler a-e scheme, i.e., an explicit scheme, when the viscosity is absent. 

Finally, an alternate space-time mesh is depicted in Fig. 30. The intrinsic geometry of 
this mesh is determined by four parameters h, r, w, and 8 with 0 < r < 1. The use of this 
mesh also results in a stencil with the shape of a tetrahedron in a 3-D space-time. In this 
Figure, mesh points marked with solid circles are centers of the SEs at the half-integer time 
levels; while those marked with open circles are centers of the SEs at the whole-integer time 
levels. Also spatial projections of the interfaces which di vide CEs are marked wi th dash 
lines. It is easy to see that each SE is associated with three CEs. According to a discussion 
given in Sec. 4 (p.46), there is an advantage that the center of each .SE be located at the 
geometric center of the top face of the union of the three CEs which are associated with 
this SE. It can be shown that the mesh has the above property if x — — 4 « 0.3589. 
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Appendix A 


The proofs for several results presented in Sec. 2 are given here. 

To prove Eqs. (2.53)-(2.58), first we shall evaluate the flux leaving each of six quadri- 
laterals that form the boundary of a CE (see Figs. 5(a) and 6(a)). As a preliminary, note 
that, in Fig. 5(a), 


area of ABGF — area of CDGB — area of EFGD = — — . (A.l) 

o 

In Fig. 6(a), we have 

area of BCG A = area of DEGC = area of FAGE = — — . (A.2) 


Equations (A.l) and (A.2) can be proved easily using the information provided in Fig. 7(a). 
Moreover, because u*(x, y, t;j , k, n ) is linear in x, y, and t (see Eq. (2.11)), its average value 
over any quadrilateral is equal to its value at the geometric center of the quadrilateral. 
With the above preparations, flux evaluation can be carried out easily using Eqs. (2.6a)- 
(2.6c), (2.9), (2.11), (A.l), and (A.2). 

For each quadrilateral, the result of flux evaluation is a formula involving a x , a y , 
(u r )y fc , and (%)”*. It can be converted to another formula involving , a+, (u+ )”*, 

and (u+ )? fc . To carry out the above conversion, note that Eqs. (2.24), (2.25), (2.29), and 
(2.30) imply that 



(A.3) 


and, for any (j, k, n ) € ft, 


/(“*)", k\ 
\ («»)",*/ 


£ 

w 


( • 

I w + b 

V h~ 


1 \ 

w — b 
~ / 


/(“<)>. A 


(A.4) 


Let (u x )”*, ( u y) 1 j t k-> ■ • ■■> be abbreviated as u x , u v , . . ., respectively. Then Eqs. (A.3) and 
(A.4) imply that 

% = ^( a J-«c)’ ( A5 ) 

ha * + (■ j - b ) a » = -jr ( a ? + 2 a i ) > ( A - 6 ) 

h ° x ~ (f + b ) ° y = it ( 2a ? + ) ’ ( A - 7 ) 
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U x 


(4.8) 

3 

y wh . 

(w -b)u+ - (w + 6)uj] , 

(A9) 

a x u x + a y u y = a+uj + a+u+, 

(4.10) 

( b 

\2 + 6j 

!«* + £«» =2uJ-u+, 

(4.11) 

/ b w\ 
\2 ~ ) 

)u x + ^u y = uj-2u+, 

(4.12) 


(A.13) 


^-d +i )“. = ^[('- 2 +^T+ ! r)“?+( fc 2 ^ 2 -T-?)<- 

■ •• - .. _ M.h; 

The conversion referred to above can be carried out using Eqs. (A.5)-(A.14). 

Consider Fig. 5(a). The results of flux evaluation involving the quadrilaterals that 
form the boundaries of CE^(j, fc, n + 1/2), t = 1, 2, 3, and (j, k,n + 1/2) 6 fii , are: 

(1) The flux leaving CE^(j, k,n + 1/2) through G'F'A'B' is 

2 wh ( , ,\ n + 1 / 2 

— (•+■?+<)*, • 

(2) The flux leaving CEj 1 ^', k,n + 1/2) through G'GFF' is 

- y|ir( a < + 2a j)[“ + 2 “< -< + +a+ti+)] 

- & K '- 2 + * - T + “H + (* + * + T - ?)<] }" +1/ • 
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(3) The flux leaving CE^(j, n + 1/2) through G'B'BG is 

T { " T ( 2a < + “’) “ _ “c + 2 < + j (“? “c + 4 “j) ] 

+ ^[( ft2+ ‘ s+ T + ^H + ( A2+ ^-T- 2 r)<]}" 

' Jl 

(4) The flux leaving CE^^j, k,n + 1/2) through AFGB is 


2wh ( + + \» 

1 u — uj — uT 

3 V < 


j+l/3,k+l/3 


(5) The flux leaving CEj 1 ^', k, n -f 1/2) through ABB' A' is 

y{^( a c + 2a J) u ~ 2u < + <-y( a c «? + <<) 

3/-* T/.o u> 2 2wb\ . /, o W 2 4wb' 

-M h+b -T + —)< + ( h+b+ T-ir, 


Hi 


i+l/3,*+l/3 


(6) The flux leaving CE^(j, k,n + 1/2) through AA’F’F is 


t{ _ t( 2 “< +a i)[ u+u t - 2u t +«)] 

^[(* , ^+T +1 rK + (*’ + *'-T- 2 rH}’ 


j+l/3,k+l/3 


(7) The flux leaving CE^O, k,n + 1/2) through G'B'C'D' is 

2wh ( _ . 4.\ n+1 / 2 

— (“- 2 “< +<)>,» • 

(8) The flux leaving CE^C?, k,n + 1/2) through G'GBB 1 is 


_^|_^( 2a + + a +) |u-u+ + 2u+ + ^(a+u+ + «)] 

+ ^[(^ +i2+ T +1 rK + ( A2+i2 -T-¥Kl} 


n+1/2 
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(9) The flux leaving CE^O'j k,n + 1/2) through G'D'DG is 

- nr { " | (°< " <) [“ - “< " < + f (“< * “c + <* “ 

N «+l/2 

+ ^ [( w + fe ) u c -( w ~ b ) u t [ 

) j,k 

(10) The flux leaving CE^O', k,n + 1/2) through CBGD is 

3 V ^ n / ;-2/3,fc+l/3 

(11) The flux leaving CE^(j, k,n + 1/2) through CDD'C' is 

+ “») u + u c' _2 “?' _ t(“<“? +a n<)] 

+ a(‘ ,+l?+ T +1 rH + (‘ ,+ ‘ , -T-¥H} 

(12) The flux leaving CE^O, k,n + 1/2) through CC'B'B is 

^{-|(°< -<) “+«< +<-x(°f u ? +a X 


+ ^ [(^ + ^)« J J ] 

J J— 2/3, 


-2/3,ifc+l/3 

(13) The flux leaving CEj 1 ^*, fc,n + 1/2) through G'D'E'F' is 

2u>/i / . _ ,\«+i/2 

— («+«< -*<)*, ■ 

(14) The flux leaving CE^^j, k, n + 1/2) through G’GDD' is 

T*{-? ( a < - <) [ u - u < - < + t(“? u 
+ ^[ (u, + 6)u ? _(u,_6)u; !’]} 


c + « 


n+l/2 

> 

i,* 


j-2/3,fc+l/3 

)] 
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(15) The flux leaving k,n + 1/2) through G'F'FG is 

+2a O u+2u <- u t + f( a < u < + a X)] 

_ ^[( A2+62_ T + ^H + ( k2+62+ T -1 r)“’]}” 

(16) The flux leaving CEg 1 ^/, k, n + 1/2) through EDGF is 

2 wh / . „ ,\ n 


n+l/2 


’ j+l/3,fc-2/3 


(17) The flux leaving CEj 1 ^,;, k,n + 1/2) through EFF'E' is 


_!f^|_|( a +_ a +) u + u + + <-^(a+u+ +aju+) 

} »» 

j'4 


j+l/3,fc-2/3 

(18) The flux leaving CE^^j, k,n + 1/2) through EE'D'D is 


~Y {^(°< + 2 <) [u-2u+ +< -^(a+u+ + «)j 


j+ 1/3, *-2/3 


Consider Fig. 6(a). The results of flux evaluation involving the quadrilaterals that 
form the boundaries of CE f(j, k, n + 1), t = 1, 2, 3, and (j, k, n + 1) 6 ^ 2 , are: 

(19) The flux leaving CE^ 2 ^/, k, n + 1) through G'C'D'E' is 

2wh ( . .\ n + 1 

~r~ u < “V if * • 

(20) The flux leaving CEj 2 ^(j, k, n + 1) through G'GCC' is 

+ 2ai 0 w ~ 2u ? + < + ^{ a t u t +«) 

} n+l 

j,k 
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util III 


(21) The flux leaving CE^^j, k, n -f 1) through G'E'EG is 


- y |~t ( 2a < + <) [“ + “<- 2 < + j- (“t “<+«)] 


(22) The flux leaving CE^\j, k, n + 1) through DCGE is 


2wh ( , \n+l/2 

— [u + uT +u+) 

3 V ' n ) j-l/3,fc-l/3 

(23) The flux leaving CE^(j, k,n + 1) through DEE'D' is 

~ + 2a 0 [ u + 2 u ?-<~x( a ? u < + °^)] 

\ n+l/2 
) j-l/3,*-l 

(24) The flux leaving CE ^\j, k, n + 1) through DD'C'C is 

+a 0 u ~ u < +2u i~T( a < u < 

^ n+l/2 

L } j- 1/3, k-1 

(25) The flux leaving CE^O, k,n + 1) through G'E'F'A' is 

2 wh ( . 

1“ (“ + 2 “< -<). t • 

(26) The flux leaving CE^\j, k, n + 1) through G'GEE 1 is 

f {'tK + <) “+“f - + j(“( n ( + *I“J) 

+S[(* 2 + 62 + T +i rK + ( A 2 +t 2 -T-T i )<]r- 


n+l/2 


n+l/2 
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(27) The flux leaving CE^O", n + 1 ) through G'A'AG is 


Hr { \ ( a ? “ a < ) f ' “ + u < + < + t ( a t‘ u < + < < ) 


} n+l 
j,k 


(28) The flux leaving CE^ (j, k, n + 1) through FEGA is 

_H|^( u _ 2u + + u +) 


n+l/2 

j+2/3,k — 1/3 


(29) The flux leaving CE^ (i, + 1) through FAA'F' is 

- t{~'t( 2g < +a 0 u -< + 2 <- j( a < u c + «) 


+ 


3/i 

tt>/i 


[ ( ^^t +4 ?H + ('* 2+42 -t-?K]} 


n+l/2 

>+2/3,*-l/3 


(30) The flux leaving CE^^j, k,n + 1) through FF'E'E is 


_u^t|^( a +_ G +) u _ u +_ u +_*t( a + u + + a + u +)j 




+ll 


n+l/2 


j+2/3,i— 1/3 


(31) The flux leaving CE 3 (j, k,n + 1 ) through G'A'B'C' is 

2 wh / . +\” +1 

I" (— <+2-J) m • 

(32) The flux leaving CEj 2 ^', k,n + 1 ) through G'GAA' is 

n+1 


} n+l 
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(33) The flux leaving CE ?\j, k,n + 1) through G'C'CG is 


- t { t(“< + 2 “» ) “ - 2u < + 4 + t ( a < u < +44)] 

' i>* 


(34) The flux leaving CE^O, k,n 4 1) through BAGC is 


2 wh / , \n+l/2 

- 2 <W 


Jfc+2/3 


(35) The flux leaving CE ?^(j,k,n 4- 1) through BCC'B' is 


ir{l(4-4) u_ 4 -4 -if (44 +44) 

'I n+1 / 2 

+ ^ [(^ + *)«< - (u> - i)u+ l 

' J-l/3,k+2/3 


(36) The flux leaving CeP(j, k,n + 1) through BB'A'A is 


f { T (°< + 2ai 0 u + 2 P ~ u t ~ j ( a c u t + «) 

-^[( ft 2 ^ 2 -T +? rH + ( ft 2 +t 2 4 -^H]}” 


-l/3,k+2/3 


With the aid of Eqs. (2.31), (2.32) and (2.35)-(2.52), Eqs. (2.53)-(2.58) are the results 
of (1)— (36) and Eqs. (2.12) and (2.13). QED. 

To prove Eqs. (2.86) and (2.87), we evaluate A^ and A^ using Eqs. (2.35)-(2.52). 
After simplifications, the results are 


a ( 1) = 3 [3(1 + i/<)( 1 + u v )(l - i/ c - i/„) + 2(1 4 u c )(l - i/ c - i/,X c 
4 2(1 4 v„)(l -v < - 4 2(1 4 ^cXl + u v)(r 

- ((cf - Uv? - Ur? + 2e C e„ + %(c(r + 2^^] , 


(A. 15) 
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and 


a (2) = 3^3(1 - ^ c )(l - v v ){\ + vc + u v ) + 2(1 - i/ c )(l + v c + i/,X c 

+ 2 (1 “ ^)(1 + I'C + + 2 (! “ I / c )( 1 ~ v n)ir ( A16 ) 

- ( 0) 2 - (£,) 2 - ( ir ? + 2 fc£, + 2 t c ( r + 20 £ r ] . 

By comparing Eqs. (2.86) and (2.87) with Eqs. (A. 15) and (A. 16), it is seen that the proof 
of the former is completed if one can prove 

-«<) 2 - «,) 2 - « r ) 2 + 2 ( ( (, + 2 £ c f r + 2 « r = (| j £) . ( 4 . 17 ) 

Proof of Eq. (A. 17): The area of the hexagon ABCDEF depicted in Fig. 7(a) is 2 wh. The 
area of ABDF depicted in the same figure is half of that. Thus 

the area of ABDF = wh. (A.18) 

Moreover, because a£, atj, and at are the lengths of the three sides of ABDF , 

the area of ABDF 

(4.19) 


Combining Eqs. (A.18)-(A.19), one has 

(aC + at] + at)(a£ + at] — at)(a£ + at — atj)(at] + at — aC) = I6w 2 h 2 . (A.20) 



A direct result of Eq. (2.32) is 

- «<) 2 - «,) 2 - Urf + 2( ( (, + 2( ( ( r + 2« r = 

[-(aC ) 4 - (at?) 4 - (at) 4 + 2(AC) 2 (Ar?) 2 + 2 (a0 2 (at) 2 + 2(a»7) 2 (at) 2 ] . 

Because 

- (AC) 4 - (A,) 4 - (at) 4 + 2(a<) 2 (a,) 2 + 2(aC) 2 (at) 2 + 2(a,) 2 (at) 2 
= (AC + AT} + At)(a£ + AT] — At)(a( + AT — Arj)(AT] + AT — A0. 


(4.21) 


(4.22) 


Equation (A. 17) is a direct result of Eqs. (A. 20) and (A. 21). QED. 
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Appendix B 


The proof for Eqs. (4.56)-(4.61) is given here. 

As a preliminary, note that Eqs. (4.22), (4.25), and (4.31) can be used to obtain 


4 

fmt ~ ~ 53 (/f<9 U 9 x ’ 


fmt = - 53 fm,t + ft,q U n) ’ 


(B. 2 ) 


In this appendix, we adopt the same convention stated following Eq. (4.37). It follows 
from Eqs. (4.34)-(4.37) that 


2 ( W ~ b W + b \ 


—h h 


m = 1,2, 3,4, 


(3.3) 


w + b w — b 


An immediate result of Eqs. (B.3) and (B.4) is 


m = 1,2, 3, 4. (B.4) 


E {fkl + fl.t “<.) = A E <( + <) • "* = 1 . 2 . 3 , 4 . (B. 5 ) 


By using Eqs. (4.18), (4.20)-(4.25), and (B.1)-(B.5), it can be shown that 




(B.6) 


( b w\ h _ . . 

~Q ) Umx 2 ~ _ 


(B.7) 


‘mi + 2 Um y ~ u inj? 


(B.8) 
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(B.9) 


** +(!-*) *-££(# 


*« -(=+*) a 

<=i 


y + -/» 
mx ' 2 J Tn V 


‘[(5 + ?)^ + 5«.]-(? + »)[(M)^ 
= SrE( 2 ^+-OK-“t). 


-f 1 + 17 * _ ^ » + ~fy 

\2 6/ Jmi 2 /my J V3 / [\2 6/ /mi 


4u>/i 

9a7 


E(/^ + 2/^)(<-2u+), 

t=l 


hfU + (I -»)/', = -^r E ('& + 2 *m) (/&<< + /£<,) • 

v y />?=! 


-Vi. + (f + ») /i. = ^3 E ( 2 /& + /&) (4:<c + /£»«) ■ 


9(a<) 


*,?=! 


fV ±—fV ±—f y 
J m g Jmi ^ Jmt 

2/i 


= 5^E (/S - /&) “< ± 4 ± «t t E ( 4H c + 4X,) 

<=i L ?=i 


and 


10 A< .. 

fl ± y/mi =F T /i, 


2/i vl 


3a< 


E (/S - /ijr) u < ± u t< ±u t ± H u tc + ft * u ™) 

<=i ?=i 


(5.10) 


(5.11) 


(5.12) 


(5.13) 


(5.14) 


(5.15) 


(5.16) 
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Note that each of Eqs. (B.15) and (B.16) represents two equations. One corresponds to 
the upper signs; while the other, to the lower signs. 

Next we shall evaluate the flux of h* m leaving each of the six quadrilaterals that form 
the boundary of a CE (see Figs. 5(a) and 6(a)). The evaluation procedure is similar to that 
described in Appendix A. For the current case, the key equations used are Eqs. (2.6a)- 
(2.6c), (4.19), (4.27)-(4.29), and (B.6)-(B.16). 

Consider Fig. 5(a). The results of flux evaluation involving the quadrilaterals that 
form the boundaries of CE^(j, k,n + 1/2), £ = 1,2, 3, and (j, k, n + 1/2) € Sl\, are: 

(1) The flux of h* m leaving CE^(j, k,n + 1/2) through G'F'A'B' is 


2wh 


(“m + U + c + U+„) 

' / Jy 


n+l/2 

k 


(2) The flux of h* m leaving CE ^\j, k,n + 1/2) through G'GFF' is 

{ E (£!< + a/I!,) [“' + 2u <! -<+E (/!>!< + /, 

W=1 L 9=1 


2 wh 
9 




n+l/2 

> 


(3) The flux of h* m leaving CE^(j, k, n + 1/2) through G' B' BG is 
2 wh 


{ E ( 2 /™!, + /»!,) f“' - + E + /?.!<) } 

l t=l 9=1 J J 

(4) The flux of h ^ leaving CE^(j, k,n + 1/2) through AFGB is 

2wh ( + + \" 

" 3 r-~ u -C- u m9j. +1/3Jk+1/3 - 

(5) The flux of leaving CE^(j, A:, n + 1/2) through ABB' A' is 


-| ^ "+1/2 
i,* 


2wh 

T~ 


|E(/m!, +2 /mT<) «< - 2»?c + «?, - E (/,.!“!<+ /<”!“!•>) ] 

U=i L 9=1 jK. 


j+l/3,k+l/3 


(6) The flux of h leaving CE^(j, k,n + 1/2) through AA'F'F is 


2 wh 
~9~ 


{ E ( 2 /?, + /*) f“« + - 2“?, - E (/,!!“!, + /, >!,) 

I <=1 L 9=1 


j+l/3,k+l/3 
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(7) The flux of h* m leaving CE^O', k, n -f 1/2) through G'B'C'D' is 

2wh ( „ 4- + \ n+1 / 2 

(8) The flux of h* m leaving CE £\j,k,n + 1/2) through G'GBB' is 


2wh 


>. n+l/2 


fe( 2 /S/ + /mt/) u < ~ u i + 2u t + E (/& u g< + ft q u m) } 

lt=i L ? =1 


(9) The flux of leaving CE £\j,k,n + 1/2) through G'D'DG is 

2wh r E (ft, - ft,) U - «tc - «i + E (ft << + } 

l 1=1 l 0=1 J ) 


?,* 

n+l/2 

► 


(10) The flux of h* m leaving CE^O, k, n + 1/2) through CBGD is 

2wh 


( [u m +2ui ( -ut v )"_ 


,*+ 1/3 

(11) The flux of leaving CE^C/, k,n + 1/2) through CDD'C' is 

-^{e (/ft, + ft,) [«. + - 2-/, - E (4X + «•*)] }” 


,-2/3,*+l/3 


(12) The flux of h ^ leaving CE^O’, k,n + 1/2) through CC'B'B is 

^{tl(ft,-ft,) «< + «£ + »?,- E (#“J< + On) } 

(13) The flux of ^ leaving CE^/, k, n + 1/2) through G'D'E'F' is 


-2/3, *+1/3 


2 wh 


(u m + «i c - 2 u+„) 


I 9 D* F * E1/ 

n+l/2 
* 


(14) The flux of leaving CE^j, Ar, n 4- 1/2) through G'GDD' is 

~ir (e (ft, - ft<) [“< - “« - “i + E (/&»£ + #“«)] 1 

l /=1 L 9=1 ) J 


n+l/2 
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(15) The flux of h* m leaving CE^^j, k,n + 1/2) through G'F'FG is 


2wh 


'4 r 4 I "4 n+1 / 2 

E (/£i+ 2 /£) “<+H-<+E (//X +/?>«) 

• <=1 L * =1 JJ y,t 


(16) The flux of h ^ leaving i,n + 1/2) through EDGF is 

2wh 


(u m — n't, + 2u+ 'j 

\ mC mn / j+l/3,fc-2/3 


(17) The flux of leaving CE^^j, fc, n + 1/2) through EFF'E' is 


2wh 


{ E (/£ - /£;) “< + + »?, - E (/£>* + OJ.) } 

l /=1 L g=l J J . 


j+l/3,*-2/3 


(18) The flux of h ^ leaving CE^\j, k,n+ 1/2) through EE'D'D is 


2wh 


i=i 


f&+ *f&) “< - K + - E (ft << + ft <) } 


j+l/3,Jfc-2/3 


Consider Fig. 6(a). The results of flux evaluation involving the quadrilaterals that 
form the boundaries of CE^(j, fc, n + 1), t — 1, 2, 3, and (j, k, n + 1) € axe: 

(19) The flux of h leaving CE^(j, &, n + 1) through G'C'D'E' is 


2wh 




n+1 

it 


(20) The flux of h* m leaving CE^(j>, k,n + 1) through G'GCC' is 
2 wh 


( 4 r 4 'i n+1 

|E(«+ 2 *&) »«-K+«S+E(^»k+^“J '0 

U=1 L q=i }) jtk 


(21) The flux of h leaving CE^(j, A:,n 1) through G'E'EG is 


2wh 


( 4 C 4 1 ) " +1 

jE( 2 /»W£) “<+<- 2 <+E(/<,X<+/<>i) } • 

U=i L g=i 


82 


(22) The flux of k* m leaving CE^(j, k, n + 1) through DCGE is 

2 wh / , _l \ ”+ 1 / 2 


( i i 

[Um+Ut f + «i,J._ 1/M _ 1/s 


(23) The flux of ^ leaving CEp^j, fc,n + 1) through DEE'D' is 
2 wh 


< 4 4 -I) "+ 1 / 2 

|E(C+<<) 

W=1 g=l J ) j-\/3,k-\/3 


(24) The flux of h* m leaving CE^\j, k, n + 1) through DD'C'C is 
2wh 


9 


{e( 2 #<+/«/) 4 + 2 <~Z (ft 4 + ft 4) } 


n+l/2 

► 

j-l/3,k-l/3 


(25) The flux of h leaving CE ^\j, k, n + 1) through G'E'F'A' is 

2wh 


(u m + 2 u+ c - u+„) 

V y Ji 


n+l 

it 


(26) The flux of leaving CE ^\j,k,n + 1) through G'GEE' is 


2wh 


f 4 r 4 1 1 " +1 

{Z( 2 ft>+ft<) »<+4~ 2 4+E(ft4+ft4) 

U=i L ?=i l) jik 


(27) The flux of h* m leaving CE ^\j,k,n + 1) through G'A'AG is 


2i oh 


E (ft> - ft<) + 4 +4,+ E (ft 4 + ft 4) 

><= i L ?=i 


n+l 




(28) The flux of leaving CE ^\j,k,n + 1) through FEGA is 


2 wh 


(29) The flux of h* m leaving CE ^\j,k,n + 1) through FAA'F' is 


(u m - 2u+ c + u+ 


n+l/2 


j+2/3,*-l/3 


nr I E ( 2 ft + ft<) I <* - «/< + 2 <, - E (ft 4 + ft 4) 


.t=i 


?=i 


n+l/2 


J ' j+2/3,fc— 1/3 
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I'll, itt 


(30) The flux of h leaving CE^O, k,n + 1) through FF'E'E is 


2 wh 


( 4 r 4 i "j n+1 / 2 

E(/ft-/ft) 

'* =1 •- ? _1 -• i+2/3,ifc-l/3 


(31) The flux of leaving CE^^j, fc,n -f 1) through G'A'B'C' is 

2 wh ( . _ . \ n+1 

(32) The flux of h* m leaving CE ^\j, k, n + 1) through G'GAA' is 

-et{ «<+<+<+£ (fti u t< + ftZ u tv) } 


n+l 


(33) The flux of h leaving CE^^j, k,n + 1) through G'C'CG is 
2wh 


( 4 r 4 ") "+ 1 / 2 

E (/ft + 2 /ft) u,-2<+u++e (/&«i+/ff*i) 

U=1 9=1 J J J,fe 


(34) The flux of leaving CE f\j, k,n + 1) through BAGC is 

2wh ( . . \ n+1 /2 

— 3-(u m + « m< -2« m ,j._ i/si+2/3 

(35) The flux of h* m leaving CE^j, k,n + 1) through BCCB' is 


2wh 

T 


r 4 r 4 ) n+1 / 2 

{ E (/ft - /Si) u < - < - «*, - E (/£“« + /£“«) f 

Ifal L 5=1 J ' j— 1/3, 

(36) The flux of leaving CE^^j, k,n + 1) through BB'A'A is 


,*+2/3 


2 wh 
~9~ 


{ E (/ft + 2 /ft) U + H - - E (/?X + /£««) 

U=i L 9=1 


n+l/2 


j-l/3,fc+2/3 


With the aid of Eqs. (4.38)-(4.55), Eqs. (4.56)-(4.61) axe the results of (l)-(36) and 
Eqs. (4.32) and (4.33). QED. 
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Appendix C 


In this appendix, first we shall prove the following theorem. 

Theorem. Let A and B be two arbitrary n x n matrices. Then AB and BA have the 
same eigenvalues. 

Proof: We have (i) 


and (ii) 


AB = A{BA)A- 


AB = B~ 1 (BA)B 


if A 1 exists; 


if B 1 exists. 


(C.l) 


(C.2) 


Thus AB ~ BA, i.e., AB is similar to BA, if either A, or B, is nonsingular. Because two 
similar matrices have the same eigenvalues [25, p.45], the proof is complete if either A, or 
B is nonsingular. 

Let both A and B be singular. Because the determinant of AB (or BA) is the product 
of the determinant of A and the determinant of B, both AB and BA axe singular. Thus 
0 is an eigenvalue for both AB and BA. We shall prove that a nonzero eigenvalue of AB 
must also be an eigenvalue of BA, and vice versa. 

Let $ be an eigenvector of AB with the eigenvalue A ^ 0. Then 


AB<j> = A <j>. 


(C.3) 


Because A ^ 0 and, by definition, $ is not a null vector, Eq. (C.3) implies that B<f> is not 
a null vector too. The last result coupled with another result of Eq. (C.3), i.e., 


BA(B<f>) = A (B<f>), 


implies that A is also an eigenvalue of BA. Conversely, it can be shown that an eigenvalue 
A ^ 0 of BA is also an eigenvalue of AB. QED. 

An immediate result of the above theorem is that the amplification matrices that 
appear on the right sides of Eqs. (5.2) and (5.3) have the same eigenvalues. Next, as a 
part of stability study, we shall investigate the two-way marching nature of the a scheme. 

According to Figs. 5(a) and 6(a), there are three CEs located immediately below any 
mesh point G‘ € ft. By definition, the mesh indices of these CEs are those of the mesh 
point G' . As shown in Figs. 12(a) and 12(b), there are also three CEs located immediately 
above any mesh point G 6 ft. However, the mesh indices of these CEs differ from those of 
the mesh point G. 

In Fig. 12(a), B' 6 ft 1? D' G fti, F' € fti, and G E ft 2 - The three CEs located imme- 
diately above point G have their mesh indices tied to points B' , D', and F' , respectively. 
In Fig. 12(b), A! € ft 2 , C' € ft 2 , E' € ft 2 , and G 6 fti- The three CEs located immediately 
above point G have their mesh indices tied to points A', C' , and E', respectively. 
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FYom Figs. 5(a), 6(a), 12(a), and 12(b), one concludes that each CE is in contact with 
(i) one and only one mesh point € fii, and (ii) one and only one mesh point G 0,2- Each 
CE is located immediately below one of these two mesh points and above the other. As an 
example, let (j, k, n + 1/2) € fti. Then (j + 1/3, k + 1/3, n) € According to Fig. 5(a), 
CE^(j, k,n + 1/2) is located immediately below the mesh point (j, fc, n + 1/2) and above 
the mesh point ( j + 1/3 , fc + 1/3, n). Note that Eq. (2.60) represents a relation among the 
marching variables associated with the above two mesh points. A similar interpretation 
can be given to each of Eqs. (2.61)-(2.65). In Sec. 2, it is shown that Eqs. (2.60)-(2.65) 
are equivalent to the defining equations of the a scheme, i.e., Eqs. (2.67)-(2.72). In the 
following, it will be shown that the former also are equivalent to the defining equations of 
the backward-marching version of the a scheme. 

Let ( j,k,n ) G f &2 and consider Fig. 12(a). Then (j — 1/3, k — 1/3, n + 1/2) G fli, 
(j + 2/3, k — 1/3, n + 1/2) G fii, and ( j — 1/3, k + 2/3, n + 1/2) G flj. By making some 
substitutions in mesh indices and exchanging expressions on the left and the right sides, 
Eqs. (2.60)-(2.62) imply that 

[«-(! + “ C 1 + v n) u t] j k = [« + (! + + C 1 + v v) u v \ /"J, fc _ 1/3 7 (C. 5) 


[ in r -i n-fl / 2 

« + (2 - v c )«f - (1 + = [« - (2 - i/()uj + (1 + ^)<j. +2/3 , (C.6) 

and 

r r , ."in+l /2 

u - (1 + i/ c )u+ + (2 - i/„X . = [« + (1 + - (2 - */,)«+ . , (C. 7) 

respectively. Let and denote the expressions on the right sides of Eqs. (C.5), 

(C.6), and (C.7), respectively. Then these equations are equivalent to 

u ”k = 5 [(! - "C - + (1 + ‘ / c)^2 1> + (1 + ‘',)4' > ] , (c.8) 

(«?)?,» = | (4 11 - 4 °) . (C.9) = 

and 

(Om = l (4" - 4") • (c. io) 

In other words, the marching variables associated with point G can be expressed in terms 
of those associated with points B\ D', and F'. . 

Let (/, fc, n + 1/2) G fti and consider Fig. 12(b). Then (j + 1/3, fc + 1/3, n + 1) G Q 25 
(j — 2/3, k + 1/3 , n + 1) G Q 2 , and (j + 1/3, & — 2/3 , n + 1) G ^ 2 - By making some 
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substitutions in mesh indices and exchanging expressions on the left and the right sides, 
Eqs. (2.63)-(2.65) imply that 


[u + (1 - */ c )uJ + (1 - */,)«+ 


n+l/2 


[« — (i — vd u + 



n+l 

i+l/3,Jt+l/3 ’ 

(C.ll) 


[u - (2 + + (1 - u n )u; 


n+l/2 


and 


u + (2 + i/{)u 


+ 

C 


~ (1 ~ v n )u+ 


■ n+l 

- >-2/3,i+l/3 ’ 

(C. 12) 


F , -i n+l/2 r 1 n+l 

[“ + (!- - ( 2 + ^XJ . 4 =[“-(!- + (2 + X“Jj /+1/3 fc _ 2/3 > 

(CU3) 

respectively. Let ij 2 \ s^\ and denote the expressions on the right sides of Eqs. (C.ll), 
(C.12), and (C.13), respectively. Then these equations are equivalent to 


n+l/2 _ 1 

3 


(i + z/ c + i/,)sS 2) + (i - ^)4 2) + (i - ^)4 2) » 


(C. 14) 


(“?)”,t 1/2 = l( J f ) -'>f ) ). (C.15) 

and 

WW-l (#>-&). (c- 16 ) 

In other words, the marching variables associated with point G can be expressed in terms 
of those associated with points A ' , C', and E'. 

The backward marching version of the a scheme consists of two marching steps. The 
first is formed by Eqs. (C.14)-(C.16); while the second is formed by Eqs. (C.8)-(C.10). To 

^(jk) 

express these steps in the forms of Eqs. (2.82) and (2.83), let the column matrices a t and 

cj(*) 

j3 t ,k = 1,2, and ^ = 1,2, 3, be defined by 


2,(1) def 1 
Or, = - 




and 




(C. 17) 


2.(1) def 1 
IT 



and 



(■£:■>) • 


(C.18) 
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C,(l) def 1 

° 3 =3 


1 + 
0 
1 


and J 1 + 1 /£ 


,<2) def 1 ( 1+ Y''” 
1 


=3 


and 


c*(2) def 1 
“2 = o 


i.(2) def 1 

“ 3 =3 


1-iV 

-1 

0 

’1-^ 

0 

-1 


and d = I -(1 - i/ c ) 


, ^ 2) def , 0 

and p 2 =1 2 + 1 /£ 


— (2 — Vff) 

1 

-C 1 ~ v n), 
1 

hi 


~J2) 


and 0 Z d = ( -(l-i/ c ) 

2 + i/, 

i.- 


(C.19) 


(C.20) 


(C.21) 


(C.22) 


Let the 3x3 matrices of rank one , fc = 1, 2, and £ = 1, 2, 3, be defined by 

(‘>¥^ ^ (k) V 


g/ ; = ^ 


(*)’ 


(C.23) 


(C.24) 


where the row matrix ( 0 t j is the transpose of the column matrix 0 t . Then Eqs. (C.8)- 
(C.10) can be expressed as 

q O’, ”) = Q^q{j - 1/3, k - 1/3, n + 1/2) + (i + 2/3, k - 1/3, n + 1/2) 

+ <> 3 1) <f(i - 1/3, +2/3,n + 1/2), (j,fe,n) G a 2 . 

Similarly, Eqs. (C.14)-(C.16) can be expressed as 

q(j , fc, n + 1/2) = g} 2) f (j + 1/3, fc + 1/3, n + 1) + Q< 2) g (j - 2/3, k + 1/3, n + 1) 

+ Q?q(j + 1/3, fc - 2/3, n + 1), (j, k, n + 1/2) € fii- 

(C.25) 

Next we shall establish several mathematical relations involving the column matrices 
<3^, 0f k \ ai \ and 0 t .To proceed, let 




(\ l + V{ l + i/, \ 

1 -(2 - l/ C ) 1 + U v 

Vl l + i/£ -(2-i/,)/ 


(C.26) 
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1 

-(1 + v ( ) 

+ 

Afg> 

1 

2 — V{ 

-(1 + u v ) 



-(1 + v c ) 

2 -u v / 


r 

-(1 - i/ c ) 



i 

2 + 

-C 1 ~ v v) 


u 

-(1 - i/ c ) 

2 + Vq / 


f 1 

i-^c 

1 - "r, \ 

Mg } 4 = 

i 

-(2 + v$) 

l-v„ 



l~ v < 

“(2 + Vr,)' 


(C. 27) 


(C.28) 


(C.29) 


Note that: (i) is the coefficient matrix of the expressions on the left sides of Eqs. (2.60) 

-(2.62). It is also that of the expressions on the right side of Eqs. (C.5)-(C.7); (ii) 
is the coefficient matrix of the expressions on the right sides of Eqs. (2.60)-(2.62). It is 
also that of the expressions on the left side of Eqs. (C.5)-(C.7); (iii) is the coefficient 
matrix of the expressions on the left sides of Eqs. (2.63)-(2.65). It is also that of the 
expressions on the right side of Eqs. (C.11)-(C.13); and (iv) is the coefficient matrix 
of the expressions on the right sides of Eqs. (2.63)-(2.65). It is also that of the expressions 
on the left side of Eqs. (C.11)-(C.13). 

Moreover, for each k = 1,2, let (i) be the 3x3 matrix formed by the column 
matrices <J^, £ = 1,2, 3; (ii) B ^ be the 3x3 matrix formed by the column matrices 0^ k \ 

£ — 1,2,3; (iii) be the 3x3 matrix formed by the column matrices a t , £ = 1,2,3; 

-Jk) 

and (iv) B (k ’ be the 3x3 matrix formed by the column matrices , £ = 1,2, 3. 

With the above preliminaries, it is easy to shown that 

M ( L k) = [^ (i) ] _1 = [i? (fc) ] * , k = 1, 2, (C.30) 

and 

M ( R k) = [i (fc) ] _1 = k = 1,2. (C.31) 

Here [A] -1 and [>!]* denote the inverse and the transpose of any matrix A. It follows from 
Eqs. (C.30) and (C.31) that 


A (k) = /, and [£(*)]* = j, fc = l,2, (C.32) 
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m 

- 


and 

[5 (t) ]‘ A (k) = I, and A (k) [F (fc) ]‘ = /, k = 1,2, (C.33) 

where I is the 3x3 identity matrix: i 

Recall the definitions of the matrices A^ k \ B^ k \ A^ k \ and B^ k \ k = 1,2. Then the 
first parts of Eqs. (C.32) and (C.33) can be expressed as 

* = 1,2, and £,€' = 1,2,3, (C.34) 


k = 1,2, and £,£' = 1,2,3, (C.35) 



respectively. Here 6tp is the Kronecker delta symbol. Also, the second parts of Eqs. (C.32) 
and (C.33) can be expressed as 



k = 1,2, 


(C.36) 


and 



* = 1 , 2 , 


(C.37) 


respectively. 

To proceed further, let 




n o 

0 ^ 




p del 

0 -1 

0 

• . . - . 

(C.38) 



\0 0 

-1/ 



We have 

P- 1 = F, 

and 

P* = F, 

(C.39) 

and 

M[ k) = A# 1 P, 

k 

= 1,2. 

(C.40) 


By taking the transpose of the expression on each side of Eq. (C.40), and using Eqs. (C.30), 
(C.31), and (C.39), one concludes that 


~JM) 


0 t = P$ fc) , * = 1,2, and £=1,2,3. 


(C.41) 
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By taking the inverse of the expression on each side of Eq. (C.40), and using Eqs. (C.30), 
(C.31), and (C.39), one concludes that 

a? = P of \ A: = 1,2, and t = 1, 2, 3. (C.42) 

Equations (C.34)-(C.37), (C.41), and (C.42) will be used in the following development. 
By using Eqs. (2.80), (C.23), (C.34), and (C.35), one has 

Q {k) Q {k) = Swaf? , * = 1,2, and = 1,2,3, (C.43) 

and 

Q {k) Q ( t k) = Sw a t k) (4 fc> ) * , k = 1,2, and tj = 1 , 2, 3. ( C.44) 

It follows from Eqs. (C.43) and (C.44) that, 

Q {k) Q {k) = Q (k) Q (k) = 0, if £ ± i'. ( k = 1, 2, and ^,£' = 1,2, 3.) (C.45) 

With the aid of Eqs. (C.36), and (C.37), Eqs. (C.43) and (C.44) also imply that 

E fffW* = E = *< k = l,2. (CM) 

e=i e=i 

Furthermore, by using Eqs. (2.80), (C.23), (C.41), and (C.42), one has 

Q {k) = PQ {k) P , and Q {k) = PQ {k) P. ( k = 1, 2, and £,£' = 1,2, 3.) (C.47) 

Equations (C.45) and (C.46) can be obtained from a more direct but less revealing 
approach. To proceed, we substitute Eq. (C.24) into the expression on the right side of 
(2.82). The result is 

S(j.*,n + 1/2) = Q^Q^qU +1, *, n + 1/2) + <?i’ ) <& 1) 9 0, * + 1, n + 1/2) 

+ QP&i'fU ~ 1, k, n + 1/2) + Q^Q^qU - 1, k + l,r> + 1/2) 

+ <?S ,) e < l 1, 9 O'. * ~ l.n + 1/2) + Q^Q^qU + 1, k - 1, n + 1/2) 

+ (<? < , 1) oS ,) + + QPWftiU,*," + 1/2). 

(C.48) 
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where (j, k,n + 1/2) € fti. By substituting Eq. (2.82) into the expression on the right side 
of Eq. (C.24), one has 


q(j,k,n) = Q^Q^qij - 1, k,n) + Q^Q^q^j, k - l,n) 

+ (j + hk,n) + Q ( 2 ) Q ( 3 1) q{j + l,k-l,n) 

+ Q^Qi^ih * + 1, n) + Q^Q^q (j — 1,& + l,n) 

+ + $? ) 02 1) + 9 O', *, n), 

where (jf, fc, n) E ^ 2 - Similarly, Eqs. (C.25) and (2.83) imply that 

9 0, k,n + l) = Q^Q^qti ~ 1, k, n + 1) + Q^Q^qU, k-l y n + l) 

+ Q^Q^qU + 1, k,n + 1) + Q% ) Q 3 ) q\ (j + 1, k - 1, n + 1) 

+ Q^Q^qQ, k + 1, n + 1) + ~ 1, k + 1 , n + 1) 

+ (Q? ) Q? ) + + Qi 2) Qi 2) )q\j,k,n + 1 ), 


(C.49) 


(C.50) 


where (jf,Ar,n + 1) € 1 ^ 2 ? and ’ , . _ 

? O', Jr, n + 1/2) = Q[ 2) Q[ n qU + 1, *, n + 1/2) + <jf ><?f 9 0, k + 1, n + 1/2) 

+ $V, 2 V(j - l,t,» + 1/2) + OfefsO' - M + l,n+ 1/2) 

+ <??’ <?f 9 0, * - 1, " + 1/2) + QfQfqii + 1, Jr - 1, n + 1/2) 

+ (Q? ) Q? ) + QfQf + qTq?) fU, M + 1/2), 

( 081 ) 

where (j, k, n+1/2) £ Oi . Because (i) and time level no can be considered as the initial time 
level, i.e., one can assign arbitrary values to the members of the set {q(j, k, n 0 ) | (j, k , n 0 ) € 
H); (ii) the column matrices in each of Eqs. (C.48)-(C.51) are associated with the same 
time level; and (iii) Eqs. (C.48)-(C.51) are valid no matter what values are assigned to the 
elements of the column matrices in these equations, a comparison of the expressions on 
the left and the right sides of each of these equations reveals that, in each of Eqs. (C.48)- 
(C.51), the first six coefficient matrices on the right side (each of them is the product 
of two matrices) must vanish; while the last coefficient matrix (the sum of three matrix 
products) must be equal to the identity matrix. As a result, Eqs. (C.45) and (C.46) follow 
from Eqs. (C.48)-(C.51). QED. 

Next Eqs. (C.45)-(C.47) will be used in a study of the amplification matrices of the 
a scheme. To proceed, let 


M (1 >(^ c ,^) d = + (C.52) 

and 

M< 2 >(0 <5 0„) = + Q< 2) e<*/ 3 >(- 2fl < + ^ + Q< 2 > e (*/ 3 )(«<-2<M (C.53) 


i 

j 
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where i = >/— 1, and — 7r < 0^,0 V < it. Equations (5.4), (5.5), (C.52), and (C.53) coupled 
with Eqs. (C.45) and (C.46) imply that 

= 1, k = 1,2. (C.54) 


In other words, M^{B^,6 V ) is the inverse of M^(B^,B V ) and vice versa. As a result, both 
axe nonsingular. Equation (C.54) can be used to show that 




Let #,,)] * and , 0^ )| , denote the component-wise complex conju- 

gates of M( k \d{,$ v ) and respectively. Then Eqs. (5.4), (5.5), (C.52), (C.53) 

coupled with Eq. (C.47) imply that 

[m<‘)(# c , #,)] = P [m<«(# c , #,)] ‘ P, k — 1,2. (C56) 

Combining Eqs. (C.39), (C.54), and (C.56), one arrives at 

m^\9 c ,9,)m^(9 c ,9,) 

Thus 

M™(9 C ,9,)M W (9 C ,0,) ~ [*<»(«<, tf,)M< 2 ) (tf c> #,)]". (C. 58) 

Let Ai , A 2 , and A 3 be the eigenvalues of the matrix on the left side of Eq. (C.58). According 
to Eq. (5.3), the last matrix is the amplification matrix for any two consecutive whole- 
integer time levels. Equation (C.58) implies that Aj, A 2 , and A 3 are also the eigenvalues 
of the matrix on the right side of this equation [25, p.45]. Because the eigenvalues of 
the matrix A*, the component- wise complex conjugate of a matrix A, are the complex 
conjugates of the eigenvalues of the matrix A , A*, A£, and A 3 , the complex conjugates of 
Ai, A 2 , and A 3 , are the eigenvalues of the matrix enclosed within the brackets on the right 
side of Eq. (C.58). According to Eq. (C.55), the last matrix is the inverse of the matrix 
on the left side of Eq. (C.58). Because the eigenvalues of the matrix A -1 , the inverse of 
a matrix A, are the reciprocals of the eigenvalues of A, one concludes that the set of A*, 
£ = 1,2,3, is identical to the set of 1 /A*, £ = 1,2,3. It does not implies that A* = 1/A*, 
£ — 1,2,3. However, it does implies that the product of AJ, £ = 1,2,3, is equal to the 
product of 1 /A*, £ = 1,2,3. As a result, we arrive at Eq. (5.7). 
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Next we shall prove that the backward characteristic projection of the mesh point 
(j, k,n + 1/2) € fti at the (n — l/2)th time level is in the interior of the n um erical domain 
of dependence of ( j , k,n + 1/2) at the same time level if and only if Eq. (2.110) is satisfied. 
In Fig. 13(c), the mesh point ( j , k,n + 1/2) is represented by point O; while its backward 
characteristic projection at the (n — l/2)th time level is represented by point P. Without 
any loss of generality, we will assume that j = k == 0. Thus 

( = fj = 0, and i = (n + 1/2) At, (C. 59) 


for point 0. Note that only the coordinates ((,ry) are given in Fig. 13(a). 

To simplify the discussion, Eq. (3.1) is converted to an equivalent form in which (, ry, 
and t are the independent variables, i.e., 


du du du 

& +a <ai + <, "a? = 0 - 


The characteristics of Eq. (C.60) is the family of straight lines defined by 


(C.60) 


£ = a{t + cj, and rj = a n t c 2 , (C.61) 

where c x and c 2 are constant along a characteristic, and vary from one characteristic to 
another. Because points O and P share the same characteristic line, Eqs. (C.59) and 
(C.61) imply that 


( = — At, rj — —a v ^t, and f = (n — l/2)At, (C.6 2) 

[ [ _ W ..1. ; . / 7 - . . : . "V , ■ ‘ ; .... • - - : 

for point P. 

The numerical domain of dependence of point O referred to in Sec. 5, a hexagon lying 
on the plane t = (n— 1/2) At, is depicted in Fig. 13(c). The coordinates ((, rj) of the vertices 
A , B, C, D, E, and F are also given in the same figure. The six edges of the hexagon and 
their equations on the (,-r) plane are 


AB 

DE 

BC 

EF 

CD 

FA 


Ary • ( + aC • rj = a(a ry, 
Ary ■ C + aC • »y = -ACAry, 
77 = A77, 


r? = -Ary, 

C = -AC, 
C = aC, 


(C.63) 


As a result, a point ((, ry) is in the interior of the hexagon ABCDEF if and only if 


|at7 • C + aC * »y| < ACAry, |ry| < A 7, and |Ci < aC- {CM) 
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Equations (C.62) and (C.64) coupled with Eqs. (2.29) and (2.31) imply that point P is in 
the interior of the hexagon ABCDEF if and only if Eq. (2.110) is satisfied. QED. 

Next we turn to the a- pi scheme. Recall that the equations used to construct the 
backward marching version of the a scheme, i.e., Eqs. (C.5)-(C.7) and (C.11)-(C.13), are 
obtained from Eqs. (2.60)-(2.65) by making substitutions in mesh indices and exchanging 
expressions on the left and the right sides. In a similar manner, the backward marching 
version of the a-pi scheme can be constructed using Eqs. (2.53)-(2.5 8) as the starting 
point. Let A^ and A^ be the determinants of the matrices formed by and 

k,i = 1,2,3, respectively. Without going into the details, it can be shown that (i) 


A (1) = 3 


3(1 + J'cXl + v,)(l u n ) - 2(1 + v ( )(l v v )U 


- 2(1 + I/,)(1 - i/ f - u n )( n - 2(1 + t^)(l + V v )£ r + | > 


and 


A< 2 > = 3 


3(1 - i/ c )(l - i/,)(l + i/c + u„) - 2(1 - i/ c )(l + u< + v n )£ c 


2 * 

- 2 (1 “ ^)(1 + "< + v T,)Zn - 2(1 - */ f )(l ~ 

and (ii) the backward marching version of the a-pi scheme exists if 

A^ ^ 0, and A^ ^ 0. 

Recall that the forward marching version exists if 

a^ / 0, and a^ ^ 0. 


(C. 65) 


(CM) 


(C.6 7) 


(C.68) 


We assume Eqs. (C.67) and (C.68). Then the matrices Q\ k \ k = 1,2, and t = 1,2,3, 
can be constructed such that the backward marching version of the a- pi scheme can also be 
represented by Eqs. (C.24) and (C.25). By using an earlier argument involving Eqs. (C.48)- 
(C.51), Eqs. (C.45) and (C.46) can also be established. Equations (C.54) and (C.55) follows 
immediately. However, it should be noted that Eq. (C.47) is not applicable for the a-fi 
scheme except for the case pi = 0. Thus Eq. (5.7) generally is not valid for the a- pi scheme. 



Appendix D 


The definition of the local CFL number v t used in Sec. 6 is given here. 

In Fig. 21, point P 0 is a mesh point (j, fc, n) (= Cl- Without any loss of generality, we 
assume that it is also the origin of the x-y plane. As explained in Sec. 5 and Appendix 
C (see Figs. 13(a)-(c)), the numerical domain of dependence of Pq on the plane with 
t = (n — 1 )a t (i.e., the (n — l)th time level) is the hexagon ABCDEF depicted in Fig. 21 . 
Because b = 0, the coordinates (x, y) of the vertices A, B, C, D, E, and F are those given 
in Fig. 21. 

The intersection of (i) the Mach cone [26, p.425] with point Pq being its vertex, and 
(ii) the plane with < = (n — l)Af, is the circle depicted in Fig. 21. Here a circle, as in the 
case of the hexagon mentioned above, implies its boundary and interior. For the Euler 
equations Eq. (4.10), and in the limit of At -4 0, this circle is the domain of dependence of 
point P 0 on the plane with t = (n — l)At. Let u, v, and c, be the x-velocity, the y- velocity, 
and the sonic speed at the point Po, respectively. Then (i) cAt is the radius of the circle, 
and (ii) x = — uAt and y = —va t axe the coordinates of the center of the circle (point Pi). 
The number i/ e will be defined such that v t < 1 if and only if the domain of dependence 
of the Euler equations (i.e., the circle) is within the interior of the numerical domain of 
dependence (i.e., the hexagon ABCDEF). 

In Fig. 21, we assume that u < 0 and v < 0. Thus x > 0 and y > 0 for point Pi . Also 
we have 

6 = f arcsin (w/y/h? +w 2 J, (0 <8 < tt/ 2) (P.l) 

and 

{ arctan |v/u|, if u ^ 0; 

7r/2, if u = 0 and v ^ 0; (^-2) 

0, if u = v = 0. 

Here we assume t hat 0 < arc tan |u/u[ < x/2. Moreover, the lengths of the line segments 
P 0 P 2 , P 0 P 3 , P 0 P 4 , and P 0 P 5 are 


IP 0 P 2 I = (c+ |u|)At, 

JP 0 P 3 J = w, 

I-P 0 P 4 I = (c+ V u 2 + v 2 cos (6 — 0 1 ) j At, 

and 

IP0P5I = 2hsin^, 

respectively. As a result of Eqs. (D.3)-(D.6), we have 

' d = l l P o P 2l _ (c + |ul)Af 

1/6 I/V5I " ^ ’ 


(P.3) 

(P.4) 

(P.5) 

(P.6) 

(P.7) 
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and 


„ def jiWl _ [c + Vu 2 + V 2 COS (6 - 0')] At 
IP0P5I 2 h sin 8 


(D. 8) 


By their definitions, (i) u e ’ < 1 if and only if the entire circle is within the domain to the 

left of the straight line A$; and (ii) v e " < 1 if and only if the entire circle is within the 

domain below the straight line B<5. Because x > 0 and y > 0 for point Po, the center of the 
circle, one concludes that the entire circle is within the interior of the hexagon ABCDEF 
if and only if v € < 1 where 

u e d = max { v t \ i/ e "}. (D.9) 


By using an argument similar to that presented above, it can be shown that, regardless 
of the signs of u and v, the entire circle is within the interior of the hexagon ABCDEF if 
and only if u e < 1 where u e is defined using Eqs. (D.1)-(D.9). 
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Table 1. — Definitions of test problems numbers 1 to 6 and the 
corresponding values of T, v ems , v em , and n c 
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Figure 1 . — The relative spatial positions of the mesh points € ft-) and the mesh 
points € fl 2 (dash lines are spatial boundaries of the conservation elements 
depicted in figs 5(a) and 5(b)). 


(-2, 4) 



Figure 2. — The spatial mesh indices (j, k) of the mesh points € fti (n = ±1/2, ±3/2, ±5/2, — ). 
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Figure 3. — The spatial mesh indices fl, k) of the mesh points 6 &2 (n = 0, ±1 , ±2, — ). 


k 



Figure 4.— The spatial mesh positions of the mesh points marked 
by • and those marked by o. 
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t 
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Figure 5. — (a) Conservation elements CEO) (j, k, n + 1/2), £ = 1, 2, 3, 
and j, k, n = 0, ±1 , ±2, •••• (b) Solution elements SEC) (j, k, n + 1/2), 
j, k, n = 0, ±1, ±2, •••). 


102 
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Figure 6. — (a) Conservation elements CE< 2 ) (j, k, n + 1), € = 1, 2, 3, 
j, k, = 1/3, 1/3 ±1,1/3 ±2, — , and n = 0, ±1 , ±2, — . (b) Solution 
elements SE( 2 ) {j, k, n + 1), j, k = 1/3, 1/3 ±1,1/3 ±2, — , and n = 0, 
± 1 . ± 2 , •••)■ 
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«4 k 4 n| 



Figure 8. — (a) The mesh points (j, k, n + 1/2), Q + 1/3, k + 1/3, n), (j - 2/3, k + 

1 /3, n) and Q + 1 /3, k - 2/3, n) with (j, k, n + 1/2) € ft-| . (b) The mesh points 
0, k, n + 1), 0 - 1/3, k - 1/3, n +1/2), (j + 2/3, k - 1/3, n + 1/2), and Q - 1/3, 
k + 2/3, n + 1/2) with 0, k, n + 1) € 112- 



Figure 9. — The stability domain of the a scheme. 
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Figure 11.— The weighted-average a-« scheme, (a) (j, k, n + 1/2) eft-), 
(b) Q, k, n + 1)€ fl 2 - 
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Rgure 14. — (a) The stability domain of a-e scheme fort = 0.1. (b) The stability domain of a-e scheme fore = 0.5. (c) The stability domain of a-e 
scheme for e = 0.8. 
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Figure 17.— {a) The stability domain of a-|i scheme for § = 10~ 5 . (b) The stability domain of a-p. scheme for| = 10 -3 . (c) The stability domain of 
a-ji scheme for | = 1(H. 


112 









iJtti 


r — Upper boundary J. 

! I 3 3 3 3 | p 


j (1,D (1,2) (1*3) (1°S) 

djl) (1,2) (1°3) (- 

1*0 • o 

(2J) (2,2) (2,3) (2, S) 


(1*S) (i,js + 1) h 


Inflow _ 
boundary 


(2, 2) (2, 3) 


(3|1) (3°2) 


(2, S) (2,jS + 1) 
o i - Outflow 


(3, 3) (3, S) 


(3, 2) (3, 3) 


(3, S) (3, S + 1) 


(R*1) (R?2) 


(R. 3) (R, S) 


(R, 2) (R, 3) 


(R, s) (R|s + 1) 


B (R + 1,1) (R + 1,2) (R + 1,3) 1 (R + 1 , S) C • • • 

O • O I • o 

(R + 1,1) (R + 1,2) (R + 1,3) i (R + 1.S) (R + 1.S + 1) 


Solid wall 


Figure 20. — The spatial locations and the new mesh indices (r, s) of mesh points (R = S = 4). 
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Figure 21 . — Definition of the CFL number. 
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E m (n), m = 1 , 2, 3, 4 




Figure 24. — Numerical results and convergence histories for 
problem #3. (a) Pressure coefficients at the mid-section of the 
computation domain (y = 0.5 in Fig. 1 9). (b) Convergence histories 
for u m , m = 1, 2, 3, 4. 
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Figure 25. — Numerical results and convergence histories for 
problem #4. (a) Pressure coefficients at the mid-section of the 
computation domain (y = 0.5 in Fig. 19). (b) Convergence historie 
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